PWGLF/FORWARD/corrections/Emperical/GenerateDispVtxCorr.C

   1 /**
   2  * Open a file
   3  *
   4  * @param filename Name of file
   5  *
   6  * @return Return value
   7  */
   8 TFile* OpenFile(const char* filename)
   9 {
  10   TFile* file = TFile::Open(filename, "READ");
  11   if (!file)
  12     Error("OpenFile", "Failed to open the file %s for reading", filename);
  13   return file;
  14 }
  15 /**
  16  * Get an object from parent list
  17  *
  18  * @param l    Parent list
  19  * @param name Name of object to find
  20  *
  21  * @return 0 in case of problems, object pointer otherwise
  22  */
  23 TObject* GetObject(const TList* l, const char* name)
  24 {
  25   if (!l) {
  26     Warning("GetObject", "No parent list given");
  27     return 0;
  28   }
  29   TObject* o = l->FindObject(name);
  30   if (!o) {
  31     Warning("GetObject", "Object named %s not found in list %s",
  32             name, l->GetName());
  33     return 0;
  34   }
  35   return o;
  36 }
  37
  38 /**
  39  * Get a list from another list
  40  *
  41  * @param l    Parent list
  42  * @param name Name of list to find
  43  *
  44  * @return 0 in case of problems, object pointer otherwise
  45  */
  46 TList* GetList(const TList* l, const char* name)
  47 {
  48   TObject* o = GetObject(l, name);
  49   if (!o) return 0;
  50   return static_cast<TList*>(o);
  51 }
  52 /**
  53  * Get a histogram from a list
  54  *
  55  * @param l    Parent list
  56  * @param name Name of histogram to find
  57  *
  58  * @return 0 in case of problems, object pointer otherwise
  59  */
  60 TH1* GetHist(const TList* l, const char* name)
  61 {
  62   TObject* o = GetObject(l, name);
  63   if (!o) return 0;
  64   return static_cast<TH1*>(o);
  65 }
  66 /**
  67  * Generate the emperical correction for unknown material in the
  68  * AliROOT geometric description.  This is done by taking the ratio
  69  * of the nominal vertex results to the results from satelitte events.
  70  *
  71  * This can be done in two ways
  72  *
  73  * - Either by using FMD results only
  74  * - Or by using the full combined FMD+V0+SPD results
  75  *
  76  * The correction is written to the file "EmpiricalCorrection.root",
  77  * which will contain a TGraphErrors for each defined centrality
  78  * class, plus a TGraphErrors object that contains the average over
  79  * all centrality classes.
  80  *
  81  * Other results should be DIVIDED by the correction obtained from
  82  * this script.
  83  *
  84  * @param fmdfmd If true, use FMD results only
  85  */
  86 void
  87 GenerateDispVtxCorr(Bool_t fmdfmd=true,
  88                     const char* fmdDispVtx="forward_dndetaDispVtx.root",
  89                     const char* fmdNomVtx="forward_dndetaNominalVtxZDC.root",
  90                     const char* fullDispVtx="dndetaPreliminaryQM12.root")
  91 {
  92
  93   TString outName    = "EmpiricalCorrection.root";
  94   TFile* fdispvtxFMD = OpenFile(fmdDispVtx);
  95   TFile* fNominalFMD = OpenFile(fmdNomVtx);
  96   TFile* fdispvtxAll = OpenFile(fullDispVtx);
  97   if (!fdispvtxFMD || !fNominalFMD || !fdispvtxAll) return;
  98
  99   TGraphErrors* corrcent[4];
 100   TMultiGraph*  mg = new TMultiGraph();
 101   mg->SetName("corrections");
 102
 103   Int_t limits[] = {     0,        5,      10,         20,    30  };
 104   Int_t colors[] = {kRed+2, kGreen+2, kBlue+1, kMagenta+1, kBlack };
 105
 106   TGraphErrors* allsym     = 0;
 107   TGraphErrors* fmddispvtx = 0;
 108   TGraphErrors* fmdnominal = 0;
 109
 110   TList* displist = static_cast<TList*>(fdispvtxFMD->Get("ForwardResults"));
 111   TList* nomlist  = static_cast<TList*>(fNominalFMD->Get("ForwardResults"));
 112
 113   Int_t nMin = 1000;
 114   // Get each defined centrality and form the ratio
 115   for(Int_t i=0; i<4; i++) {
 116     Int_t cl = limits[i];
 117     Int_t ch = limits[i+1];
 118     corrcent[i] = new TGraphErrors;
 119     corrcent[i]->SetName(Form("correction_cent_%03d_%03d",cl, ch));
 120     corrcent[i]->SetTitle(Form("%2d%% - %2d%%", cl, ch));
 121     corrcent[i]->GetHistogram()->SetXTitle("#eta");
 122     corrcent[i]->GetHistogram()->SetYTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
 123                                            "dN_{ch}/d#eta|_{satelitte}}");
 124     corrcent[i]->SetMarkerColor(colors[i]);
 125     corrcent[i]->SetLineColor(colors[i]);
 126     corrcent[i]->SetMarkerStyle(8);
 127     corrcent[i]->SetFillStyle(0);
 128     corrcent[i]->SetFillColor(0);
 129
 130     mg->Add(corrcent[i]);
 131     TGraphErrors* allsym =
 132       static_cast<TGraphErrors*>(fdispvtxAll->Get(Form("graphErrors_cent_%d_%d",
 133                                                        cl, ch)));
 134
 135     TString folderName   = Form("cent%03d_%03d",cl, ch);
 136     TList*  dispcentlist = GetList(displist, folderName);
 137     TList*  nomcentlist  = GetList(nomlist, folderName);
 138     TString histName     = "dndetaForward_rebin05";
 139     TH1D*   hDisp        = GetHist(dispcentlist, histName);
 140     TH1D*   hNominal     = GetHist(nomcentlist, histName);
 141
 142     // Make our temporary graph
 143     if   (fmdfmd) fmddispvtx = new TGraphErrors(hDisp);
 144     else          fmddispvtx = static_cast<TGraphErrors*>(allsym->Clone());
 145     fmdnominal               = new TGraphErrors(hNominal);
 146
 147     Int_t nPoints = 0;
 148
 149     for(Int_t n=0; n<fmdnominal->GetN(); n++) {
 150       Double_t eta        = fmdnominal->GetX()[n];
 151       Double_t nommult    = fmdnominal->GetY()[n];
 152       Double_t nommulterr = fmdnominal->GetErrorY(n);
 153       Double_t etaerr     = fmdnominal->GetErrorX(n);
 154
 155       // Ignore empty bins
 156       if(nommult < 0.0001) continue;
 157
 158       // Find the corresponding bin from the dispaclaced vertex analysis
 159       for(Int_t m=0; m < fmddispvtx->GetN(); m++) {
 160         Double_t eta1        = fmddispvtx->GetX()[m];
 161         Double_t dispmult    = fmddispvtx->GetY()[m];
 162         Double_t dispmulterr = fmddispvtx->GetErrorY(m);
 163
 164         if(TMath::Abs(eta-eta1) >= 0.001) continue;
 165
 166         Double_t corr  = nommult/dispmult;
 167         Double_t rd    = dispmulterr/dispmult;
 168         Double_t rn    = nommulterr/nommult;
 169         Double_t error = (1/corr)*TMath::Sqrt(TMath::Power(rd,2) +
 170                                               TMath::Power(rn,2));
 171         corrcent[i]->SetPoint(nPoints,eta,corr);
 172         corrcent[i]->SetPointError(nPoints,etaerr,error);
 173         nPoints++;
 174
 175       }
 176       //std::cout<<eta<<"  "<<nommult<<std::endl;
 177     }
 178     nMin = TMath::Min(nPoints, nMin);
 179   }
 180
 181   // Calulate the average
 182   TGraphErrors* average = new TGraphErrors;
 183   average->SetName(Form("correction_average"));
 184   average->SetTitle(Form("%2d%% - %2d%%", limits[0], limits[4]));
 185   average->SetMarkerColor(colors[4]);
 186   average->SetLineColor(colors[4]);
 187   average->SetMarkerStyle(8);
 188   average->SetFillStyle(0);
 189   average->SetFillColor(0);
 190   mg->Add(average);
 191
 192   for(Int_t k = 0; k < nMin; k++) {
 193     Double_t mean   = 0;
 194     Double_t error2 = 0;
 195     Double_t sumw2  = 0;
 196
 197     // Loop over centralities
 198     for(Int_t l=0; l<4; l++) {
 199       Double_t eta  = corrcent[l]->GetX()[k];
 200       Double_t corr = corrcent[l]->GetY()[k];
 201       Double_t err  = corrcent[l]->GetErrorY(k);
 202       Double_t err2 = err * err;
 203       Double_t w    = 1 / err2;
 204       sumw2         += w;
 205       mean          += w * corr;
 206 #if 0
 207       mean   += corr;
 208       error2 += TMath::Power(err,2);
 209 #endif
 210     }
 211     mean           /= sumw2;
 212     Double_t error = TMath::Sqrt(1. / sumw2);
 213 #if 0
 214     mean           /= 4;
 215     Double_t error =  TMath::Sqrt(error2) / 4;
 216 #endif
 217
 218     average->SetPoint(k,eta,mean);
 219     average->SetPointError(k,0.125,error);
 220   }
 221
 222   Double_t min = +1000;
 223   Double_t max = -1000;
 224   TMultiGraph* ratios = new TMultiGraph;
 225   for(Int_t l=0; l<4; l++) {
 226     Int_t cl        = limits[l];
 227     Int_t ch        = limits[l+1];
 228     TGraphErrors* r =
 229       static_cast<TGraphErrors*>(corrcent[l]->Clone(Form("ratio%02d%02d",
 230                                                         cl,ch)));
 231     ratios->Add(r);
 232     for (Int_t k = 0; k < r->GetN(); k++) {
 233       Double_t x = r->GetX()[k];
 234       Double_t y = r->GetY()[k];
 235       Double_t a = average->Eval(x);
 236       Double_t s = (y-a)/a;
 237       r->SetPoint(k, x, s);
 238       min  = TMath::Min(s, min);
 239       max  = TMath::Max(s, max);
 240     }
 241   }
 242   Printf("Min=%f max=%f", min, max);
 243   // Draw the results
 244   TCanvas* canvas = new TCanvas("overview","overview",800,600);
 245   TPad* overview =  new TPad("overview", "Overview", 0, .3, 1, 1, 0, 0);
 246   overview->SetBottomMargin(0);
 247   overview->SetTopMargin(0.02);
 248   overview->SetRightMargin(0.02);
 249   overview->Draw();
 250   overview->SetGridx();
 251   overview->cd();
 252
 253   mg->Draw("APEL");
 254   mg->GetXaxis()->SetTitle("#eta");
 255   mg->GetYaxis()->SetTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
 256                            "dN_{ch}/d#eta|_{satelitte}}");
 257   overview->Clear();
 258   mg->DrawClone("APEC");
 259   TLegend* leg = overview->BuildLegend(0.35,0.6,0.6,0.975, "Centrality");
 260   leg->SetFillColor(0);
 261   leg->SetFillStyle(0);
 262   leg->SetBorderSize(0);
 263
 264   canvas->cd();
 265   TPad* details =  new TPad("details", "Details", 0, 0, 1, .3, 0, 0);
 266   details->SetBottomMargin(0.15);
 267   details->SetTopMargin(0);
 268   details->SetRightMargin(0.02);
 269   details->Draw();
 270   details->SetGridx();
 271   details->cd();
 272
 273
 274   ratios->Draw("APEL");
 275   ratios->GetXaxis()->SetTitle("#eta");
 276   ratios->GetYaxis()->SetTitle("#frac{c-#LTc#GT}{#LTc#GT}");
 277   ratios->GetYaxis()->SetTitleSize(0.09);
 278   ratios->GetYaxis()->SetTitleOffset(0.5);
 279   ratios->GetYaxis()->SetLabelSize(0.08);
 280   ratios->GetXaxis()->SetTitleSize(0.09);
 281   ratios->GetXaxis()->SetTitleOffset(0.5);
 282   ratios->GetXaxis()->SetLabelSize(0.08);
 283   ratios->GetYaxis()->SetNdivisions(10);
 284
 285   Double_t x1 = ratios->GetXaxis()->GetXmin();
 286   Double_t x2 = ratios->GetXaxis()->GetXmax();
 287   Double_t y1 = min; // ratios->GetHistogram()->GetMinimum();
 288   Double_t y2 = max; // ratios->GetHistogram()->GetMaximum();
 289   TGraphAsymmErrors* band = new TGraphAsymmErrors(2);
 290   band->SetName("band");
 291   band->SetTitle(Form("Error band min=%4.2f%% max=%4.2f%%", -100*min, 100*max));
 292   band->SetPoint(0, x1, 0);
 293   band->SetPoint(1, x2, 0);
 294   band->SetPointError(0, 0, 0, -y1, y2);
 295   band->SetPointError(1, 0, 0, -y1, y2);
 296   band->SetFillColor(kYellow+2);
 297   band->SetFillStyle(3001);
 298   band->SetLineStyle(2);
 299   band->SetLineColor(kBlack);
 300   band->Draw("A E3 L");
 301   band->GetXaxis()->SetTitle("#eta");
 302   band->GetYaxis()->SetTitle("#frac{c-#LTc#GT}{#LTc#GT}");
 303   band->GetYaxis()->SetTitleSize(0.09);
 304   band->GetYaxis()->SetTitleOffset(0.5);
 305   band->GetYaxis()->SetLabelSize(0.08);
 306   band->GetXaxis()->SetTitleSize(0.09);
 307   band->GetXaxis()->SetTitleOffset(0.5);
 308   band->GetXaxis()->SetLabelSize(0.08);
 309   band->GetYaxis()->SetNdivisions(10);
 310
 311   details->Clear();
 312   // band->DrawClone("A E3 L");
 313   ratios->GetListOfGraphs()->AddFirst(band, "E3 L");
 314   ratios->DrawClone("PEC same");
 315
 316   leg = details->BuildLegend(0.35,0.6,0.6,0.975, "Centrality");
 317   leg->SetFillColor(0);
 318   leg->SetFillStyle(0);
 319   leg->SetBorderSize(0);
 320
 321   // Store the results
 322   TFile* fout = TFile::Open(outName,"RECREATE");
 323
 324   for(Int_t p=0; p<4; p++) corrcent[p]->Write();
 325   average->Write("average");
 326   mg->Write("all");
 327   ratios->Write("all");
 328   fout->ls();
 329   fout->Close();
 330   Info("GenerateDispVtxCorr", "Corrections written to %s", outName.Data());
 331 }
 332 //
 333 // EOF
 334 //