| | 1 | #include <TH1.h> |
| | 2 | #include <THStack.h> |
| | 3 | #include <TGraphErrors.h> |
| | 4 | #include <TGraphAsymmErrors.h> |
| | 5 | #include <TMultiGraph.h> |
| | 6 | #include <TFile.h> |
| | 7 | #include <TList.h> |
| | 8 | #include <TString.h> |
| | 9 | #include <TError.h> |
| | 10 | #include <TSystem.h> |
| | 11 | #include <TROOT.h> |
| | 12 | #include <TMath.h> |
| | 13 | #include <TCanvas.h> |
| | 14 | #include <TPad.h> |
| | 15 | #include <TStyle.h> |
| | 16 | #include <TLegend.h> |
| | 17 | #include <TLegendEntry.h> |
| | 18 | #include <TLatex.h> |
| | 19 | #include <TImage.h> |
| | 20 | |
| | 21 | /** |
| | 22 | * Class to draw dN/deta results |
| | 23 | * |
| | 24 | */ |
| | 25 | struct dNdetaDrawer |
| | 26 | { |
| | 27 | /** |
| | 28 | * POD of data for range zooming |
| | 29 | */ |
| | 30 | struct RangeParam |
| | 31 | { |
| | 32 | TAxis* fMasterAxis; // Master axis |
| | 33 | TAxis* fSlave1Axis; // First slave axis |
| | 34 | TVirtualPad* fSlave1Pad; // First slave pad |
| | 35 | TAxis* fSlave2Axis; // Second slave axis |
| | 36 | TVirtualPad* fSlave2Pad; // Second slave pad |
| | 37 | }; |
| | 38 | //__________________________________________________________________ |
| | 39 | /** |
| | 40 | * Constructor |
| | 41 | * |
| | 42 | */ |
| | 43 | dNdetaDrawer() |
| | 44 | : fShowOthers(false), // Bool_t |
| | 45 | fShowAlice(false), // Bool_t |
| | 46 | fShowRatios(false), // Bool_t |
| | 47 | fShowLeftRight(false), // Bool_t |
| | 48 | fRebin(5), // UShort_t |
| | 49 | fCutEdges(false), // Bool_t |
| | 50 | fTitle(""), // TString |
| | 51 | fHHDFile(""), // TString |
| | 52 | fTrigString(0), // TNamed* |
| | 53 | fSNNString(0), // TNamed* |
| | 54 | fSysString(0), // TNamed* |
| | 55 | fVtxAxis(0), // TAxis* |
| | 56 | fForward(0), // TH1* |
| | 57 | fForwardMC(0), // TH1* |
| | 58 | fForwardHHD(0), // TH1* |
| | 59 | fTruth(0), // TH1* |
| | 60 | fCentral(0), // TH1* |
| | 61 | fForwardSym(0), // TH1* |
| | 62 | fForwardMCSym(0), // TH1* |
| | 63 | fForwardHHDSym(0), // TH1* |
| | 64 | fTriggers(0), // TH1* |
| | 65 | fRangeParam(0) |
| | 66 | |
| | 67 | { |
| | 68 | fRangeParam = new RangeParam; |
| | 69 | fRangeParam->fMasterAxis = 0; |
| | 70 | fRangeParam->fSlave1Axis = 0; |
| | 71 | fRangeParam->fSlave1Pad = 0; |
| | 72 | fRangeParam->fSlave2Axis = 0; |
| | 73 | fRangeParam->fSlave2Pad = 0; |
| | 74 | } |
| | 75 | //================================================================== |
| | 76 | /** |
| | 77 | * @{ |
| | 78 | * @name Set parameters |
| | 79 | */ |
| | 80 | /** |
| | 81 | * Show other (UA5, CMS, ...) data |
| | 82 | * |
| | 83 | * @param x Whether to show or not |
| | 84 | */ |
| | 85 | void SetShowOthers(Bool_t x) { fShowOthers = x; } |
| | 86 | //__________________________________________________________________ |
| | 87 | /** |
| | 88 | * Show ALICE published data |
| | 89 | * |
| | 90 | * @param x Wheter to show or not |
| | 91 | */ |
| | 92 | void SetShowAlice(Bool_t x) { fShowAlice = x; } |
| | 93 | //__________________________________________________________________ |
| | 94 | /** |
| | 95 | * Whether to show ratios or not. If there's nothing to compare to, |
| | 96 | * the ratio panel will be implicitly disabled |
| | 97 | * |
| | 98 | * @param x Whether to show or not |
| | 99 | */ |
| | 100 | void SetShowRatios(Bool_t x) { fShowRatios = x; } |
| | 101 | //__________________________________________________________________ |
| | 102 | /** |
| | 103 | * |
| | 104 | * Whether to show the left/right asymmetry |
| | 105 | * |
| | 106 | * @param x To show or not |
| | 107 | */ |
| | 108 | void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; } |
| | 109 | //__________________________________________________________________ |
| | 110 | /** |
| | 111 | * Set the rebinning factor |
| | 112 | * |
| | 113 | * @param x Rebinning factor (must be a divisor in the number of bins) |
| | 114 | */ |
| | 115 | void SetRebin(UShort_t x) { fRebin = x; } |
| | 116 | //__________________________________________________________________ |
| | 117 | /** |
| | 118 | * Wheter to cut away the edges |
| | 119 | * |
| | 120 | * @param x Whether or not to cut away edges |
| | 121 | */ |
| | 122 | void SetCutEdges(Bool_t x) { fCutEdges = x; } |
| | 123 | //__________________________________________________________________ |
| | 124 | /** |
| | 125 | * Set the title of the plot |
| | 126 | * |
| | 127 | * @param x Title |
| | 128 | */ |
| | 129 | void SetTitle(TString x) { fTitle = x; } |
| | 130 | //__________________________________________________________________ |
| | 131 | /** |
| | 132 | * Set the file name of the file containing the HHD results |
| | 133 | * |
| | 134 | * @param fn File name |
| | 135 | */ |
| | 136 | void SetHHDFile(const char* fn) { fHHDFile = fn; } |
| | 137 | /* @} */ |
| | 138 | //================================================================== |
| | 139 | /** |
| | 140 | * @{ |
| | 141 | * @name Override settings from input |
| | 142 | */ |
| | 143 | /** |
| | 144 | * Override setting from file |
| | 145 | * |
| | 146 | * @param sNN Center of mass energy per nucleon pair (GeV) |
| | 147 | */ |
| | 148 | void SetSNN(UShort_t sNN) |
| | 149 | { |
| | 150 | fSNNSTring = new TNamed("sNN", Form("%04dGeV", sNN)); |
| | 151 | fSNNString->SetUniqueID(sNN); |
| | 152 | } |
| | 153 | //__________________________________________________________________ |
| | 154 | /** |
| | 155 | * Set the collision system |
| | 156 | * - 1: pp |
| | 157 | * - 2: PbPb |
| | 158 | * |
| | 159 | * @param sys collision system |
| | 160 | */ |
| | 161 | void SetSys(UShort_t sys) |
| | 162 | { |
| | 163 | fSNNString = new TNamed("sys", (sys == 1 ? "pp" : |
| | 164 | sys == 2 ? "PbPb" : "unknown")); |
| | 165 | fSNNString->SetUniqueID(sys); |
| | 166 | } |
| | 167 | //__________________________________________________________________ |
| | 168 | /** |
| | 169 | * Set the vertex range in centimeters |
| | 170 | * |
| | 171 | * @param vzMin Min @f$ v_z@f$ |
| | 172 | * @param vzMax Max @f$ v_z@f$ |
| | 173 | */ |
| | 174 | void SetVertexRange(Double_t vzMin, Double_t vzMax) |
| | 175 | { |
| | 176 | fVtxAxis = new TAxis(10, vzMin, vzMax); |
| | 177 | fVtxAxis->SetName("vtxAxis"); |
| | 178 | fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax)); |
| | 179 | } |
| | 180 | //__________________________________________________________________ |
| | 181 | void SetTrigger(UShort_t trig) |
| | 182 | { |
| | 183 | fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" : |
| | 184 | trig & 0x2 ? "INEL>0" : |
| | 185 | trig & 0x4 ? "NSD" : |
| | 186 | "unknown")); |
| | 187 | fTrigString->SetUniqueID(trig); |
| | 188 | } |
| | 189 | |
| | 190 | |
| | 191 | //================================================================== |
| | 192 | /** |
| | 193 | * @{ |
| | 194 | * @name Main steering functions |
| | 195 | */ |
| | 196 | /** |
| | 197 | * Run the code to produce the final result. |
| | 198 | * |
| | 199 | * @param filename File containing the data |
| | 200 | */ |
| | 201 | void Run(const char* filename="forward_dndeta.C") |
| | 202 | { |
| | 203 | if (!Open(filename)) return; |
| | 204 | |
| | 205 | Double_t max = 0; |
| | 206 | |
| | 207 | // Create our stack of results |
| | 208 | THStack* results = StackResults(max); |
| | 209 | |
| | 210 | // Create our stack of other results |
| | 211 | TMultiGraph* other = 0; |
| | 212 | if (fShowOthers || fShowAlice) other = StackOther(max); |
| | 213 | |
| | 214 | Double_t smax = 0; |
| | 215 | THStack* ratios = 0; |
| | 216 | if (fShowRatios) ratios = StackRatios(other, smax); |
| | 217 | |
| | 218 | Double_t amax = 0; |
| | 219 | THStack* leftright = 0; |
| | 220 | if (fShowLeftRight) leftright = StackLeftRight(amax); |
| | 221 | |
| | 222 | Plot(results, other, max, ratios, smax, leftright, amax); |
| | 223 | } |
| | 224 | |
| | 225 | //__________________________________________________________________ |
| | 226 | /** |
| | 227 | * Open input file, and find data |
| | 228 | * |
| | 229 | * @param filename File name |
| | 230 | * |
| | 231 | * @return true on success |
| | 232 | */ |
| | 233 | Bool_t Open(const char* filename) |
| | 234 | { |
| | 235 | TFile* file = TFile::Open(filename, "READ"); |
| | 236 | if (!file) { |
| | 237 | Error("Open", "Cannot open %s", filename); |
| | 238 | return false; |
| | 239 | } |
| | 240 | |
| | 241 | TList* results = static_cast<TList*>(file->Get("ForwardResults")); |
| | 242 | if (!results) { |
| | 243 | Error("Open", "Couldn't find list ForwardResults"); |
| | 244 | return false; |
| | 245 | } |
| | 246 | |
| | 247 | fForward = GetResult(results, "dndetaForward"); |
| | 248 | fForwardMC = GetResult(results, "dndetaForwardMC"); |
| | 249 | fTruth = GetResult(results, "dndetaTruth"); |
| | 250 | fCentral = GetResult(results, "dndetaCentral"); |
| | 251 | |
| | 252 | if (!fTrigString) |
| | 253 | fTrigString = static_cast<TNamed*>(results->FindObject("trigString")); |
| | 254 | if (!fSNNString) |
| | 255 | fSNNString = static_cast<TNamed*>(results->FindObject("sNN")); |
| | 256 | if (!fSysString) |
| | 257 | fSysString = static_cast<TNamed*>(results->FindObject("sys")); |
| | 258 | if (!fVtxAxis) |
| | 259 | fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis")); |
| | 260 | |
| | 261 | if (!fTrigString) fTrigString = new TNamed("trigString", "unknown"); |
| | 262 | if (!fSNNString) fSNNString = new TNamed("sNN", "unknown"); |
| | 263 | if (!fSysString) fSysString = new TNamed("sys", "unknown"); |
| | 264 | if (!fVtxAxis) { |
| | 265 | fVtxAxis = new TAxis(1,0,0); |
| | 266 | fVtxAxis->SetName("vtxAxis"); |
| | 267 | fVtxAxis->SetTitle("v_{z} range unspecified"); |
| | 268 | } |
| | 269 | |
| | 270 | TList* sums = static_cast<TList*>(file->Get("ForwardSums")); |
| | 271 | if (sums) |
| | 272 | fTriggers = GetResult(sums, "triggers"); |
| | 273 | |
| | 274 | if (!fForward) { |
| | 275 | Error("Open", "Couldn't find the result of the forward analysis"); |
| | 276 | return false; |
| | 277 | } |
| | 278 | file->Close(); |
| | 279 | |
| | 280 | |
| | 281 | fForwardHHD = GetHHD(); |
| | 282 | |
| | 283 | return true; |
| | 284 | } |
| | 285 | //__________________________________________________________________ |
| | 286 | /** |
| | 287 | * Make a histogram stack of results |
| | 288 | * |
| | 289 | * @param max On return, the maximum value in the stack |
| | 290 | * |
| | 291 | * @return Newly allocated stack |
| | 292 | */ |
| | 293 | THStack* StackResults(Double_t& max) |
| | 294 | { |
| | 295 | THStack* stack = new THStack("results", "Stack of Results"); |
| | 296 | max = TMath::Max(max, AddHistogram(stack, fTruth, "e5 p")); |
| | 297 | max = TMath::Max(max, AddHistogram(stack, fForwardHHD, "", fForwardHHDSym)); |
| | 298 | max = TMath::Max(max, AddHistogram(stack, fForwardMC, "", fForwardMCSym)); |
| | 299 | max = TMath::Max(max, AddHistogram(stack, fCentral, "")); |
| | 300 | max = TMath::Max(max, AddHistogram(stack, fForward, "", fForwardSym)); |
| | 301 | return stack; |
| | 302 | } |
| | 303 | //__________________________________________________________________ |
| | 304 | /** |
| | 305 | * Make a histogram stack of results |
| | 306 | * |
| | 307 | * @param max On return, the maximum value in the stack |
| | 308 | * |
| | 309 | * @return Newly allocated stack |
| | 310 | */ |
| | 311 | TMultiGraph* StackOther(Double_t& max) const |
| | 312 | { |
| | 313 | gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C"); |
| | 314 | Int_t error = 0; |
| | 315 | Bool_t onlya = (fShowOthers ? false : true); |
| | 316 | Int_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0); |
| | 317 | UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0); |
| | 318 | Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d);", |
| | 319 | snn,trg,onlya)); |
| | 320 | if (error) { |
| | 321 | Error("StackOther", "Failed to execute GetData(%d,%d,%d)", |
| | 322 | snn, trg, onlya); |
| | 323 | return 0; |
| | 324 | } |
| | 325 | if (!ret) { |
| | 326 | Warning("StackOther", "No other data found for sNN=%d, trigger=%d", |
| | 327 | snn, trg); |
| | 328 | return 0; |
| | 329 | } |
| | 330 | TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret); |
| | 331 | |
| | 332 | TGraphAsymmErrors* o = 0; |
| | 333 | TIter next(other->GetListOfGraphs()); |
| | 334 | while ((o = static_cast<TGraphAsymmErrors*>(next()))) |
| | 335 | max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY())); |
| | 336 | |
| | 337 | return other; |
| | 338 | } |
| | 339 | //__________________________________________________________________ |
| | 340 | /** |
| | 341 | * Make a histogram stack of ratios of results to other data |
| | 342 | * |
| | 343 | * @param max On return, the maximum value in the stack |
| | 344 | * |
| | 345 | * @return Newly allocated stack |
| | 346 | */ |
| | 347 | THStack* StackRatios(TMultiGraph* others, Double_t& max) |
| | 348 | { |
| | 349 | THStack* ratios = new THStack("ratios", "Ratios"); |
| | 350 | |
| | 351 | if (others) { |
| | 352 | TGraphAsymmErrors* ua5_1 = 0; |
| | 353 | TGraphAsymmErrors* ua5_2 = 0; |
| | 354 | TGraphAsymmErrors* alice = 0; |
| | 355 | TGraphAsymmErrors* cms = 0; |
| | 356 | TGraphAsymmErrors* o = 0; |
| | 357 | TIter nextG(others->GetListOfGraphs()); |
| | 358 | while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) { |
| | 359 | ratios->Add(Ratio(fForward, o, max)); |
| | 360 | ratios->Add(Ratio(fForwardSym, o, max)); |
| | 361 | ratios->Add(Ratio(fForwardHHD, o, max)); |
| | 362 | ratios->Add(Ratio(fForwardHHDSym, o, max)); |
| | 363 | ratios->Add(Ratio(fCentral, o, max)); |
| | 364 | TString oName(o->GetName()); |
| | 365 | oName.ToLower(); |
| | 366 | if (oName.Contains("ua5")) { if (ua5_1) ua5_2 = o; else ua5_1 = o; } |
| | 367 | if (oName.Contains("alice")) alice = o; |
| | 368 | if (oName.Contains("cms")) cms = o; |
| | 369 | } |
| | 370 | if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max)); |
| | 371 | if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max)); |
| | 372 | if (cms && alice) ratios->Add(Ratio(alice, cms, max)); |
| | 373 | } |
| | 374 | |
| | 375 | // If we have data from HHD's analysis, then do the ratio of |
| | 376 | // our result to that data. |
| | 377 | if (fForwardHHD) { |
| | 378 | ratios->Add(Ratio(fForward, fForwardHHD, max)); |
| | 379 | ratios->Add(Ratio(fForwardSym, fForwardHHDSym, max)); |
| | 380 | } |
| | 381 | |
| | 382 | // Do comparison to MC |
| | 383 | if (fForwardMC) { |
| | 384 | ratios->Add(Ratio(fForward, fForwardMC, max)); |
| | 385 | ratios->Add(Ratio(fForwardSym, fForwardMCSym, max)); |
| | 386 | } |
| | 387 | |
| | 388 | // Check if we have ratios |
| | 389 | if (!ratios->GetHists() || |
| | 390 | (ratios->GetHists()->GetEntries() <= 0)) { |
| | 391 | delete ratios; |
| | 392 | ratios = 0; |
| | 393 | } |
| | 394 | return ratios; |
| | 395 | } |
| | 396 | //__________________________________________________________________ |
| | 397 | /** |
| | 398 | * Make a histogram stack of the left-right asymmetry |
| | 399 | * |
| | 400 | * @param max On return, the maximum value in the stack |
| | 401 | * |
| | 402 | * @return Newly allocated stack |
| | 403 | */ |
| | 404 | THStack* StackLeftRight(Double_t& max) |
| | 405 | { |
| | 406 | THStack* ret = new THStack("leftright", "Left-right asymmetry"); |
| | 407 | ret->Add(Asymmetry(fForward, max)); |
| | 408 | ret->Add(Asymmetry(fForwardHHD, max)); |
| | 409 | ret->Add(Asymmetry(fForwardMC, max)); |
| | 410 | |
| | 411 | if (!ret->GetHists() || |
| | 412 | (ret->GetHists()->GetEntries() <= 0)) { |
| | 413 | delete ret; |
| | 414 | ret = 0; |
| | 415 | } |
| | 416 | return ret; |
| | 417 | } |
| | 418 | //__________________________________________________________________ |
| | 419 | /** |
| | 420 | * Plot the results |
| | 421 | * |
| | 422 | * @param results Results |
| | 423 | * @param others Other data |
| | 424 | * @param max Max value |
| | 425 | * @param ratios Stack of ratios (optional) |
| | 426 | * @param rmax Maximum diviation from 1 of ratios |
| | 427 | * @param leftright Stack of left-right asymmetry (optional) |
| | 428 | * @param amax Maximum diviation from 1 of asymmetries |
| | 429 | */ |
| | 430 | void Plot(THStack* results, |
| | 431 | TMultiGraph* others, |
| | 432 | Double_t max, |
| | 433 | THStack* ratios, |
| | 434 | Double_t rmax, |
| | 435 | THStack* leftright, |
| | 436 | Double_t amax) |
| | 437 | { |
| | 438 | gStyle->SetOptTitle(0); |
| | 439 | gStyle->SetTitleFont(132, "xyz"); |
| | 440 | gStyle->SetLabelFont(132, "xyz"); |
| | 441 | |
| | 442 | Int_t h = 800; |
| | 443 | Int_t w = 800; // h / TMath::Sqrt(2); |
| | 444 | if (!ratios) w *= 1.4; |
| | 445 | if (!leftright) w *= 1.4; |
| | 446 | TCanvas* c = new TCanvas("Results", "Results", w, h); |
| | 447 | c->SetFillColor(0); |
| | 448 | c->SetBorderSize(0); |
| | 449 | c->SetBorderMode(0); |
| | 450 | |
| | 451 | Double_t y1 = 0; |
| | 452 | Double_t y2 = 0; |
| | 453 | Double_t y3 = 0; |
| | 454 | if (ratios) y1 = 0.3; |
| | 455 | if (leftright) { |
| | 456 | if (y1 > 0.0001) { |
| | 457 | y2 = 0.2; |
| | 458 | y1 = 0.4; |
| | 459 | } |
| | 460 | else { |
| | 461 | y1 = 0.2; |
| | 462 | y2 = y1; |
| | 463 | } |
| | 464 | } |
| | 465 | PlotResults(results, others, max, y1); |
| | 466 | c->cd(); |
| | 467 | |
| | 468 | PlotRatios(ratios, rmax, y2, y1); |
| | 469 | c->cd( ); |
| | 470 | |
| | 471 | PlotLeftRight(leftright, amax, y3, y2); |
| | 472 | c->cd(); |
| | 473 | |
| | 474 | |
| | 475 | Int_t vMin = fVtxAxis->GetXmin(); |
| | 476 | Int_t vMax = fVtxAxis->GetXmax(); |
| | 477 | TString trg(fTrigString->GetTitle()); |
| | 478 | Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX()); |
| | 479 | trg = trg.Strip(TString::kBoth); |
| | 480 | TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev", |
| | 481 | fSysString->GetTitle(), |
| | 482 | fSNNString->GetTitle(), |
| | 483 | trg.Data(), |
| | 484 | vMin < 0 ? 'm' : 'p', TMath::Abs(vMin), |
| | 485 | vMax < 0 ? 'm' : 'p', TMath::Abs(vMax), |
| | 486 | nev)); |
| | 487 | c->SaveAs(Form("%s.png", base.Data())); |
| | 488 | c->SaveAs(Form("%s.root", base.Data())); |
| | 489 | c->SaveAs(Form("%s.C", base.Data())); |
| | 490 | } |
| | 491 | //__________________________________________________________________ |
| | 492 | /** |
| | 493 | * Plot the results |
| | 494 | * |
| | 495 | * @param results Results |
| | 496 | * @param others Other data |
| | 497 | * @param max Maximum |
| | 498 | * @param yd Bottom position of pad |
| | 499 | */ |
| | 500 | void PlotResults(THStack* results, TMultiGraph* others, |
| | 501 | Double_t max, Double_t yd) |
| | 502 | { |
| | 503 | // Make a sub-pad for the result itself |
| | 504 | TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0); |
| | 505 | p1->SetTopMargin(0.05); |
| | 506 | p1->SetBorderSize(0); |
| | 507 | p1->SetBorderMode(0); |
| | 508 | p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1); |
| | 509 | p1->SetRightMargin(0.05); |
| | 510 | p1->SetGridx(); |
| | 511 | p1->SetTicks(1,1); |
| | 512 | p1->SetNumber(1); |
| | 513 | p1->Draw(); |
| | 514 | p1->cd(); |
| | 515 | |
| | 516 | results->SetMaximum(1.15*max); |
| | 517 | results->SetMinimum(yd > 0.00001 ? -0.1 : 0); |
| | 518 | |
| | 519 | FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}"); |
| | 520 | |
| | 521 | p1->Clear(); |
| | 522 | results->DrawClone("nostack e1"); |
| | 523 | |
| | 524 | fRangeParam->fSlave1Axis = results->GetXaxis(); |
| | 525 | fRangeParam->fSlave1Pad = p1; |
| | 526 | |
| | 527 | // Draw other data |
| | 528 | if (others) { |
| | 529 | TGraphAsymmErrors* o = 0; |
| | 530 | TIter nextG(others->GetListOfGraphs()); |
| | 531 | while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) |
| | 532 | o->DrawClone("same p"); |
| | 533 | } |
| | 534 | |
| | 535 | // Make a legend in the main result pad |
| | 536 | TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35); |
| | 537 | l->SetNColumns(2); |
| | 538 | l->SetFillColor(0); |
| | 539 | l->SetFillStyle(0); |
| | 540 | l->SetBorderSize(0); |
| | 541 | l->SetTextFont(132); |
| | 542 | |
| | 543 | // Put a title on top |
| | 544 | TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data()); |
| | 545 | tit->SetNDC(); |
| | 546 | tit->SetTextFont(132); |
| | 547 | tit->SetTextSize(0.05); |
| | 548 | tit->Draw(); |
| | 549 | |
| | 550 | // Put a nice label in the plot |
| | 551 | TString eS; |
| | 552 | UShort_t snn = fSNNString->GetUniqueID(); |
| | 553 | const char* sys = fSysString->GetTitle(); |
| | 554 | if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000); |
| | 555 | else eS = Form("%3dGeV", snn); |
| | 556 | TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", |
| | 557 | sys, |
| | 558 | eS.Data(), |
| | 559 | fTrigString->GetTitle())); |
| | 560 | tt->SetNDC(); |
| | 561 | tt->SetTextFont(132); |
| | 562 | tt->SetTextAlign(33); |
| | 563 | tt->Draw(); |
| | 564 | |
| | 565 | // Put number of accepted events on the plot |
| | 566 | Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX()); |
| | 567 | TLatex* et = new TLatex(.93, .83, Form("%d events", nev)); |
| | 568 | et->SetNDC(); |
| | 569 | et->SetTextFont(132); |
| | 570 | et->SetTextAlign(33); |
| | 571 | et->Draw(); |
| | 572 | |
| | 573 | // Put number of accepted events on the plot |
| | 574 | if (fVtxAxis) { |
| | 575 | TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle()); |
| | 576 | vt->SetNDC(); |
| | 577 | vt->SetTextFont(132); |
| | 578 | vt->SetTextAlign(33); |
| | 579 | vt->Draw(); |
| | 580 | } |
| | 581 | // results->Draw("nostack e1 same"); |
| | 582 | |
| | 583 | fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName()); |
| | 584 | fRangeParam->fSlave1Pad = p1; |
| | 585 | |
| | 586 | |
| | 587 | // Mark the plot as preliminary |
| | 588 | TLatex* pt = new TLatex(.12, .93, "Preliminary"); |
| | 589 | pt->SetNDC(); |
| | 590 | pt->SetTextFont(22); |
| | 591 | pt->SetTextSize(0.07); |
| | 592 | pt->SetTextColor(kRed+1); |
| | 593 | pt->SetTextAlign(13); |
| | 594 | pt->Draw(); |
| | 595 | |
| | 596 | if (!gSystem->AccessPathName("ALICE.png")) { |
| | 597 | TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0); |
| | 598 | logo->SetFillStyle(0); |
| | 599 | logo->Draw(); |
| | 600 | logo->cd(); |
| | 601 | TImage* i = TImage::Create(); |
| | 602 | i->ReadImage("ALICE.png"); |
| | 603 | i->Draw(); |
| | 604 | } |
| | 605 | p1->cd(); |
| | 606 | } |
| | 607 | //__________________________________________________________________ |
| | 608 | /** |
| | 609 | * Plot the ratios |
| | 610 | * |
| | 611 | * @param ratios Ratios to plot (if any) |
| | 612 | * @param max Maximum diviation from 1 |
| | 613 | * @param y1 Lower y coordinate of pad |
| | 614 | * @param y2 Upper y coordinate of pad |
| | 615 | */ |
| | 616 | void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2) |
| | 617 | { |
| | 618 | if (!ratios) return; |
| | 619 | bool isBottom = (y1 < 0.0001); |
| | 620 | Double_t yd = y2 - y1; |
| | 621 | // Make a sub-pad for the result itself |
| | 622 | TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0); |
| | 623 | p2->SetTopMargin(0.001); |
| | 624 | p2->SetRightMargin(0.05); |
| | 625 | p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001); |
| | 626 | p2->SetGridx(); |
| | 627 | p2->SetTicks(1,1); |
| | 628 | p2->SetNumber(2); |
| | 629 | p2->Draw(); |
| | 630 | p2->cd(); |
| | 631 | |
| | 632 | // Fix up axis |
| | 633 | FixAxis(ratios, 1/yd/1.7, "Ratios", 7); |
| | 634 | |
| | 635 | ratios->SetMaximum(1+TMath::Max(.22,1.05*max)); |
| | 636 | ratios->SetMinimum(1-TMath::Max(.32,1.05*max)); |
| | 637 | p2->Clear(); |
| | 638 | ratios->DrawClone("nostack e1"); |
| | 639 | |
| | 640 | |
| | 641 | // Make a legend |
| | 642 | TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9, |
| | 643 | isBottom ? .6 : .4); |
| | 644 | l2->SetNColumns(2); |
| | 645 | l2->SetFillColor(0); |
| | 646 | l2->SetFillStyle(0); |
| | 647 | l2->SetBorderSize(0); |
| | 648 | l2->SetTextFont(132); |
| | 649 | |
| | 650 | // Make a nice band from 0.9 to 1.1 |
| | 651 | TGraphErrors* band = new TGraphErrors(2); |
| | 652 | band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1); |
| | 653 | band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1); |
| | 654 | band->SetPointError(0, 0, .1); |
| | 655 | band->SetPointError(1, 0, .1); |
| | 656 | band->SetFillColor(kYellow+2); |
| | 657 | band->SetFillStyle(3002); |
| | 658 | band->SetLineStyle(2); |
| | 659 | band->SetLineWidth(1); |
| | 660 | band->Draw("3 same"); |
| | 661 | band->DrawClone("X L same"); |
| | 662 | |
| | 663 | // Replot the ratios on top |
| | 664 | ratios->DrawClone("nostack e1 same"); |
| | 665 | |
| | 666 | if (isBottom) { |
| | 667 | fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName()); |
| | 668 | p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", |
| | 669 | fRangeParam)); |
| | 670 | } |
| | 671 | else { |
| | 672 | fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName()); |
| | 673 | fRangeParam->fSlave2Pad = p2; |
| | 674 | } |
| | 675 | } |
| | 676 | //__________________________________________________________________ |
| | 677 | /** |
| | 678 | * Plot the asymmetries |
| | 679 | * |
| | 680 | * @param ratios Asymmetries to plot (if any) |
| | 681 | * @param max Maximum diviation from 1 |
| | 682 | * @param y1 Lower y coordinate of pad |
| | 683 | * @param y2 Upper y coordinate of pad |
| | 684 | */ |
| | 685 | void PlotLeftRight(THStack* leftright, Double_t max, |
| | 686 | Double_t y1, Double_t y2) |
| | 687 | { |
| | 688 | if (!leftright) return; |
| | 689 | bool isBottom = (y1 < 0.0001); |
| | 690 | Double_t yd = y2 - y1; |
| | 691 | // Make a sub-pad for the result itself |
| | 692 | TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0); |
| | 693 | p3->SetTopMargin(0.001); |
| | 694 | p3->SetRightMargin(0.05); |
| | 695 | p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001); |
| | 696 | p3->SetGridx(); |
| | 697 | p3->SetTicks(1,1); |
| | 698 | p3->SetNumber(2); |
| | 699 | p3->Draw(); |
| | 700 | p3->cd(); |
| | 701 | |
| | 702 | TH1* dummy = 0; |
| | 703 | if (leftright->GetHists()->GetEntries() == 1) { |
| | 704 | // Add dummy histogram |
| | 705 | dummy = new TH1F("dummy","", 10, -6, 6); |
| | 706 | dummy->SetLineColor(0); |
| | 707 | dummy->SetFillColor(0); |
| | 708 | dummy->SetMarkerColor(0); |
| | 709 | leftright->Add(dummy); |
| | 710 | } |
| | 711 | |
| | 712 | // Fix up axis |
| | 713 | FixAxis(leftright, 1/yd/1.7, "Right/Left", 4); |
| | 714 | |
| | 715 | leftright->SetMaximum(1+TMath::Max(.12,1.05*max)); |
| | 716 | leftright->SetMinimum(1-TMath::Max(.15,1.05*max)); |
| | 717 | p3->Clear(); |
| | 718 | leftright->DrawClone("nostack e1"); |
| | 719 | |
| | 720 | |
| | 721 | // Make a legend |
| | 722 | TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5); |
| | 723 | l2->SetNColumns(2); |
| | 724 | l2->SetFillColor(0); |
| | 725 | l2->SetFillStyle(0); |
| | 726 | l2->SetBorderSize(0); |
| | 727 | l2->SetTextFont(132); |
| | 728 | #ifndef __CINT__ |
| | 729 | if (dummy) { |
| | 730 | TList* prims = l2->GetListOfPrimitives(); |
| | 731 | TIter next(prims); |
| | 732 | TLegendEntry* o = 0; |
| | 733 | while ((o = static_cast<TLegendEntry*>(next()))) { |
| | 734 | TString lbl(o->GetLabel()); |
| | 735 | if (lbl != "dummy") continue; |
| | 736 | prims->Remove(o); |
| | 737 | break; |
| | 738 | } |
| | 739 | } |
| | 740 | #endif |
| | 741 | // Make a nice band from 0.9 to 1.1 |
| | 742 | TGraphErrors* band = new TGraphErrors(2); |
| | 743 | band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1); |
| | 744 | band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1); |
| | 745 | band->SetPointError(0, 0, .05); |
| | 746 | band->SetPointError(1, 0, .05); |
| | 747 | band->SetFillColor(kYellow+2); |
| | 748 | band->SetFillStyle(3002); |
| | 749 | band->SetLineStyle(2); |
| | 750 | band->SetLineWidth(1); |
| | 751 | band->Draw("3 same"); |
| | 752 | band->DrawClone("X L same"); |
| | 753 | |
| | 754 | leftright->DrawClone("nostack e1 same"); |
| | 755 | if (isBottom) { |
| | 756 | fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName()); |
| | 757 | p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", |
| | 758 | fRangeParam)); |
| | 759 | } |
| | 760 | else { |
| | 761 | fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName()); |
| | 762 | fRangeParam->fSlave2Pad = p3; |
| | 763 | } |
| | 764 | } |
| | 765 | /** @} */ |
| | 766 | //================================================================== |
| | 767 | /** |
| | 768 | * @{ |
| | 769 | * @name Data utility functions |
| | 770 | */ |
| | 771 | /** |
| | 772 | * Get a result from the passed list |
| | 773 | * |
| | 774 | * @param list List to search |
| | 775 | * @param name Object name to search for |
| | 776 | * |
| | 777 | * @return |
| | 778 | */ |
| | 779 | TH1* GetResult(TList* list, const char* name) const |
| | 780 | { |
| | 781 | if (!list) return 0; |
| | 782 | |
| | 783 | TH1* ret = static_cast<TH1*>(list->FindObject(name)); |
| | 784 | if (!ret) |
| | 785 | Warning("GetResult", "Histogram %s not found", name); |
| | 786 | |
| | 787 | return ret; |
| | 788 | } |
| | 789 | //__________________________________________________________________ |
| | 790 | /** |
| | 791 | * Get the result from previous analysis code |
| | 792 | * |
| | 793 | * @param fn File to open |
| | 794 | * @param nsd Whether this is NSD |
| | 795 | * |
| | 796 | * @return null or result of previous analysis code |
| | 797 | */ |
| | 798 | TH1* GetHHD() |
| | 799 | { |
| | 800 | if (fHHDFile.IsNull()) return 0; |
| | 801 | const char* fn = fHHDFile.Data(); |
| | 802 | Bool_t nsd = (fTrigString ? fTrigString->GetUniqueID() & 0x4 : false); |
| | 803 | TDirectory* savdir = gDirectory; |
| | 804 | if (gSystem->AccessPathName(fn)) { |
| | 805 | Warning("GetHHD", "Output of HHD analysis (%s) not available", fn); |
| | 806 | return 0; |
| | 807 | } |
| | 808 | TFile* file = TFile::Open(fn, "READ"); |
| | 809 | if (!file) { |
| | 810 | Warning("GetHHD", "couldn't open HHD file %s", fn); |
| | 811 | return 0; |
| | 812 | } |
| | 813 | TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : ""))); |
| | 814 | TH1* h = static_cast<TH1*>(file->Get(hist.Data())); |
| | 815 | if (!h) { |
| | 816 | Warning("GetHHD", "Couldn't find HHD histogram %s in %s", |
| | 817 | hist.Data(), fn); |
| | 818 | file->Close(); |
| | 819 | savdir->cd(); |
| | 820 | return 0; |
| | 821 | } |
| | 822 | TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd")); |
| | 823 | r->SetTitle("ALICE Forward (HHD)"); |
| | 824 | r->SetFillStyle(0); |
| | 825 | r->SetFillColor(0); |
| | 826 | r->SetMarkerStyle(21); |
| | 827 | r->SetMarkerColor(kPink+1); |
| | 828 | r->SetDirectory(0); |
| | 829 | |
| | 830 | file->Close(); |
| | 831 | savdir->cd(); |
| | 832 | return r; |
| | 833 | } |
| | 834 | //__________________________________________________________________ |
| | 835 | /** |
| | 836 | * Add a histogram to the stack after possibly rebinning it |
| | 837 | * |
| | 838 | * @param stack Stack to add to |
| | 839 | * @param hist histogram |
| | 840 | * @param option Draw options |
| | 841 | * |
| | 842 | * @return Maximum of histogram |
| | 843 | */ |
| | 844 | Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const |
| | 845 | { |
| | 846 | // Check if we have input |
| | 847 | if (!hist) return 0; |
| | 848 | |
| | 849 | // Rebin if needed |
| | 850 | Rebin(hist); |
| | 851 | |
| | 852 | stack->Add(hist, option); |
| | 853 | return hist->GetMaximum(); |
| | 854 | } |
| | 855 | //__________________________________________________________________ |
| | 856 | /** |
| | 857 | * Add a histogram to the stack after possibly rebinning it |
| | 858 | * |
| | 859 | * @param stack Stack to add to |
| | 860 | * @param hist histogram |
| | 861 | * @param option Draw options |
| | 862 | * @param sym On return, the data symmetriced (added to stack) |
| | 863 | * |
| | 864 | * @return Maximum of histogram |
| | 865 | */ |
| | 866 | Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option, |
| | 867 | TH1*& sym) const |
| | 868 | { |
| | 869 | // Check if we have input |
| | 870 | if (!hist) return 0; |
| | 871 | |
| | 872 | // Rebin if needed |
| | 873 | Rebin(hist); |
| | 874 | stack->Add(hist, option); |
| | 875 | |
| | 876 | // Now symmetrice the histogram |
| | 877 | sym = Symmetrice(hist); |
| | 878 | stack->Add(sym, option); |
| | 879 | |
| | 880 | return hist->GetMaximum(); |
| | 881 | } |
| | 882 | |
| | 883 | //__________________________________________________________________ |
| | 884 | /** |
| | 885 | * Rebin a histogram |
| | 886 | * |
| | 887 | * @param h Histogram to rebin |
| | 888 | * @param rebin Rebinning factor |
| | 889 | * |
| | 890 | * @return |
| | 891 | */ |
| | 892 | virtual void Rebin(TH1* h) const |
| | 893 | { |
| | 894 | if (fRebin <= 1) return; |
| | 895 | |
| | 896 | Int_t nBins = h->GetNbinsX(); |
| | 897 | if(nBins % fRebin != 0) { |
| | 898 | Warning("Rebin", "Rebin factor %d is not a devisor of current number " |
| | 899 | "of bins %d in the histogram %s", fRebin, nBins, h->GetName()); |
| | 900 | return; |
| | 901 | } |
| | 902 | |
| | 903 | // Make a copy |
| | 904 | TH1* tmp = static_cast<TH1*>(h->Clone("tmp")); |
| | 905 | tmp->Rebin(fRebin); |
| | 906 | tmp->SetDirectory(0); |
| | 907 | |
| | 908 | // The new number of bins |
| | 909 | Int_t nBinsNew = nBins / fRebin; |
| | 910 | for(Int_t i = 1;i<= nBinsNew; i++) { |
| | 911 | Double_t content = 0; |
| | 912 | Double_t sumw = 0; |
| | 913 | Double_t wsum = 0; |
| | 914 | Int_t nbins = 0; |
| | 915 | for(Int_t j = 1; j<=fRebin;j++) { |
| | 916 | Int_t bin = (i-1)*fRebin + j; |
| | 917 | Double_t c = h->GetBinContent(bin); |
| | 918 | |
| | 919 | if (c <= 0) continue; |
| | 920 | |
| | 921 | if (fCutEdges) { |
| | 922 | if (h->GetBinContent(bin+1)<=0 || |
| | 923 | h->GetBinContent(bin-1)) { |
| | 924 | Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", |
| | 925 | bin, c, h->GetName(), |
| | 926 | bin+1, h->GetBinContent(bin+1), |
| | 927 | bin-1, h->GetBinContent(bin-1)); |
| | 928 | continue; |
| | 929 | } |
| | 930 | } |
| | 931 | Double_t e = h->GetBinError(bin); |
| | 932 | Double_t w = 1 / (e*e); // 1/c/c |
| | 933 | content += c; |
| | 934 | sumw += w; |
| | 935 | wsum += w * c; |
| | 936 | nbins++; |
| | 937 | } |
| | 938 | |
| | 939 | if(content > 0 && nbins > 0) { |
| | 940 | tmp->SetBinContent(i, wsum / sumw); |
| | 941 | tmp->SetBinError(i,1./TMath::Sqrt(sumw)); |
| | 942 | } |
| | 943 | } |
| | 944 | |
| | 945 | // Finally, rebin the histogram, and set new content |
| | 946 | h->Rebin(fRebin); |
| | 947 | h->Reset(); |
| | 948 | for(Int_t i = 1; i<= nBinsNew; i++) { |
| | 949 | h->SetBinContent(i,tmp->GetBinContent(i)); |
| | 950 | h->SetBinError(i, tmp->GetBinError(i)); |
| | 951 | } |
| | 952 | |
| | 953 | delete tmp; |
| | 954 | } |
| | 955 | //__________________________________________________________________ |
| | 956 | /** |
| | 957 | * Make an extension of @a h to make it symmetric about 0 |
| | 958 | * |
| | 959 | * @param h Histogram to symmertrice |
| | 960 | * |
| | 961 | * @return Symmetric extension of @a h |
| | 962 | */ |
| | 963 | TH1* Symmetrice(const TH1* h) const |
| | 964 | { |
| | 965 | Int_t nBins = h->GetNbinsX(); |
| | 966 | TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName()))); |
| | 967 | s->SetTitle(Form("%s (mirrored)", h->GetTitle())); |
| | 968 | s->Reset(); |
| | 969 | s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin()); |
| | 970 | s->SetMarkerColor(h->GetMarkerColor()); |
| | 971 | s->SetMarkerSize(h->GetMarkerSize()); |
| | 972 | s->SetMarkerStyle(h->GetMarkerStyle()+4); |
| | 973 | s->SetFillColor(h->GetFillColor()); |
| | 974 | s->SetFillStyle(h->GetFillStyle()); |
| | 975 | s->SetDirectory(0); |
| | 976 | |
| | 977 | // Find the first and last bin with data |
| | 978 | Int_t first = nBins+1; |
| | 979 | Int_t last = 0; |
| | 980 | for (Int_t i = 1; i <= nBins; i++) { |
| | 981 | if (h->GetBinContent(i) <= 0) continue; |
| | 982 | first = TMath::Min(first, i); |
| | 983 | last = TMath::Max(last, i); |
| | 984 | } |
| | 985 | |
| | 986 | Double_t xfirst = h->GetBinCenter(first-1); |
| | 987 | Int_t f1 = h->GetXaxis()->FindBin(-xfirst); |
| | 988 | Int_t l2 = s->GetXaxis()->FindBin(xfirst); |
| | 989 | for (Int_t i = f1, j=l2; i <= last; i++,j--) { |
| | 990 | s->SetBinContent(j, h->GetBinContent(i)); |
| | 991 | s->SetBinError(j, h->GetBinError(i)); |
| | 992 | } |
| | 993 | // Fill in overlap bin |
| | 994 | s->SetBinContent(l2+1, h->GetBinContent(first)); |
| | 995 | s->SetBinError(l2+1, h->GetBinError(first)); |
| | 996 | return s; |
| | 997 | } |
| | 998 | //__________________________________________________________________ |
| | 999 | /** |
| | 1000 | * Calculate the left-right asymmetry of input histogram |
| | 1001 | * |
| | 1002 | * @param h Input histogram |
| | 1003 | * @param max On return, the maximum distance from 1 of the histogram |
| | 1004 | * |
| | 1005 | * @return Asymmetry |
| | 1006 | */ |
| | 1007 | TH1* Asymmetry(TH1* h, Double_t& max) |
| | 1008 | { |
| | 1009 | if (!h) return 0; |
| | 1010 | |
| | 1011 | TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName()))); |
| | 1012 | // Int_t oBins = h->GetNbinsX(); |
| | 1013 | // Double_t high = h->GetXaxis()->GetXmax(); |
| | 1014 | // Double_t low = h->GetXaxis()->GetXmin(); |
| | 1015 | // Double_t dBin = (high - low) / oBins; |
| | 1016 | // Int_t tBins = Int_t(2*high/dBin+.5); |
| | 1017 | // ret->SetBins(tBins, -high, high); |
| | 1018 | ret->Reset(); |
| | 1019 | ret->SetTitle(Form("%s (+/-)", h->GetTitle())); |
| | 1020 | ret->SetYTitle("Right/Left"); |
| | 1021 | Int_t nBins = h->GetNbinsX(); |
| | 1022 | for (Int_t i = 1; i <= nBins; i++) { |
| | 1023 | Double_t x = h->GetBinCenter(i); |
| | 1024 | if (x > 0) break; |
| | 1025 | |
| | 1026 | Double_t c1 = h->GetBinContent(i); |
| | 1027 | Double_t e1 = h->GetBinError(i); |
| | 1028 | if (c1 <= 0) continue; |
| | 1029 | |
| | 1030 | Int_t j = h->FindBin(-x); |
| | 1031 | if (j <= 0 || j > nBins) continue; |
| | 1032 | |
| | 1033 | Double_t c2 = h->GetBinContent(j); |
| | 1034 | Double_t e2 = h->GetBinError(j); |
| | 1035 | |
| | 1036 | Double_t c12 = c1*c1; |
| | 1037 | Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12)); |
| | 1038 | |
| | 1039 | Int_t k = ret->FindBin(x); |
| | 1040 | ret->SetBinContent(k, c2/c1); |
| | 1041 | ret->SetBinError(k, e); |
| | 1042 | } |
| | 1043 | max = TMath::Max(max, RatioMax(ret)); |
| | 1044 | |
| | 1045 | return ret; |
| | 1046 | } |
| | 1047 | //__________________________________________________________________ |
| | 1048 | /** |
| | 1049 | * Transform a graph into a histogram |
| | 1050 | * |
| | 1051 | * @param g |
| | 1052 | * |
| | 1053 | * @return |
| | 1054 | */ |
| | 1055 | TH1* Graph2Hist(const TGraphAsymmErrors* g) const |
| | 1056 | { |
| | 1057 | Int_t nBins = g->GetN(); |
| | 1058 | TArrayF bins(nBins+1); |
| | 1059 | Double_t dx = 0; |
| | 1060 | for (Int_t i = 0; i < nBins; i++) { |
| | 1061 | Double_t x = g->GetX()[i]; |
| | 1062 | Double_t exl = g->GetEXlow()[i]; |
| | 1063 | Double_t exh = g->GetEXhigh()[i]; |
| | 1064 | bins.fArray[i] = x-exl; |
| | 1065 | bins.fArray[i+1] = x+exh; |
| | 1066 | Double_t dxi = exh+exl; |
| | 1067 | if (i == 0) dx = dxi; |
| | 1068 | else if (dxi != dx) dx = 0; |
| | 1069 | } |
| | 1070 | TString name(g->GetName()); |
| | 1071 | TString title(g->GetTitle()); |
| | 1072 | TH1D* h = 0; |
| | 1073 | if (dx != 0) { |
| | 1074 | h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]); |
| | 1075 | } |
| | 1076 | else { |
| | 1077 | h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray); |
| | 1078 | } |
| | 1079 | h->SetMarkerStyle(g->GetMarkerStyle()); |
| | 1080 | h->SetMarkerColor(g->GetMarkerColor()); |
| | 1081 | h->SetMarkerSize(g->GetMarkerSize()); |
| | 1082 | |
| | 1083 | return h; |
| | 1084 | } |
| | 1085 | /* @} */ |
| | 1086 | //================================================================== |
| | 1087 | /** |
| | 1088 | * @{ |
| | 1089 | * @name Ratio utility functions |
| | 1090 | */ |
| | 1091 | /** |
| | 1092 | * Get the maximum diviation from 1 in the passed ratio |
| | 1093 | * |
| | 1094 | * @param h Ratio histogram |
| | 1095 | * |
| | 1096 | * @return Max diviation from 1 |
| | 1097 | */ |
| | 1098 | Double_t RatioMax(TH1* h) const |
| | 1099 | { |
| | 1100 | Int_t nBins = h->GetNbinsX(); |
| | 1101 | Double_t ret = 0; |
| | 1102 | for (Int_t i = 1; i <= nBins; i++) { |
| | 1103 | Double_t c = h->GetBinContent(i); |
| | 1104 | if (c == 0) continue; |
| | 1105 | Double_t e = h->GetBinError(i); |
| | 1106 | Double_t d = TMath::Abs(1-c-e); |
| | 1107 | ret = TMath::Max(d, ret); |
| | 1108 | } |
| | 1109 | return ret; |
| | 1110 | } |
| | 1111 | //__________________________________________________________________ |
| | 1112 | /** |
| | 1113 | * Compute the ratio of @a h to @a g. @a g is evaluated at the bin |
| | 1114 | * centers of @a h |
| | 1115 | * |
| | 1116 | * @param h Numerator |
| | 1117 | * @param g Divisor |
| | 1118 | * |
| | 1119 | * @return h/g |
| | 1120 | */ |
| | 1121 | TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const |
| | 1122 | { |
| | 1123 | if (!h || !g) return 0; |
| | 1124 | |
| | 1125 | TH1* ret = static_cast<TH1*>(h->Clone("tmp")); |
| | 1126 | ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName())); |
| | 1127 | ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle())); |
| | 1128 | ret->Reset(); |
| | 1129 | ret->SetMarkerStyle(g->GetMarkerStyle()); |
| | 1130 | ret->SetMarkerColor(h->GetMarkerColor()); |
| | 1131 | ret->SetMarkerSize(0.9*g->GetMarkerSize()); |
| | 1132 | Double_t xlow = g->GetX()[0]; |
| | 1133 | Double_t xhigh = g->GetX()[g->GetN()-1]; |
| | 1134 | if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; } |
| | 1135 | |
| | 1136 | for (Int_t i = 1; i <= h->GetNbinsX(); i++) { |
| | 1137 | Double_t c = h->GetBinContent(i); |
| | 1138 | if (c <= 0) continue; |
| | 1139 | |
| | 1140 | Double_t x = h->GetBinCenter(i); |
| | 1141 | if (x < xlow || x > xhigh) continue; |
| | 1142 | |
| | 1143 | Double_t f = g->Eval(x); |
| | 1144 | if (f <= 0) continue; |
| | 1145 | |
| | 1146 | ret->SetBinContent(i, c / f); |
| | 1147 | ret->SetBinError(i, h->GetBinError(i) / f); |
| | 1148 | } |
| | 1149 | if (ret->GetEntries() <= 0) { |
| | 1150 | delete ret; |
| | 1151 | ret = 0; |
| | 1152 | } |
| | 1153 | else |
| | 1154 | max = TMath::Max(RatioMax(ret), max); |
| | 1155 | return ret; |
| | 1156 | } |
| | 1157 | //__________________________________________________________________ |
| | 1158 | /** |
| | 1159 | * Make the ratio of h1 to h2 |
| | 1160 | * |
| | 1161 | * @param h1 First histogram (numerator) |
| | 1162 | * @param h2 Second histogram (denominator) |
| | 1163 | * |
| | 1164 | * @return h1 / h2 |
| | 1165 | */ |
| | 1166 | TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const |
| | 1167 | { |
| | 1168 | if (!h1 || !h2) return 0; |
| | 1169 | TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s", |
| | 1170 | h1->GetName(), |
| | 1171 | h2->GetName()))); |
| | 1172 | t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle())); |
| | 1173 | t1->Divide(h2); |
| | 1174 | t1->SetMarkerColor(h1->GetMarkerColor()); |
| | 1175 | t1->SetMarkerStyle(h2->GetMarkerStyle()); |
| | 1176 | max = TMath::Max(RatioMax(t1), max); |
| | 1177 | return t1; |
| | 1178 | } |
| | 1179 | //__________________________________________________________________ |
| | 1180 | /** |
| | 1181 | * Calculate the ratio of two graphs - g1 / g2 |
| | 1182 | * |
| | 1183 | * @param g1 Numerator |
| | 1184 | * @param g2 Denominator |
| | 1185 | * |
| | 1186 | * @return g1 / g2 in a histogram |
| | 1187 | */ |
| | 1188 | TH1* Ratio(const TGraphAsymmErrors* g1, |
| | 1189 | const TGraphAsymmErrors* g2, Double_t& max) const |
| | 1190 | { |
| | 1191 | Int_t nBins = g1->GetN(); |
| | 1192 | TArrayF bins(nBins+1); |
| | 1193 | Double_t dx = 0; |
| | 1194 | for (Int_t i = 0; i < nBins; i++) { |
| | 1195 | Double_t x = g1->GetX()[i]; |
| | 1196 | Double_t exl = g1->GetEXlow()[i]; |
| | 1197 | Double_t exh = g1->GetEXhigh()[i]; |
| | 1198 | bins.fArray[i] = x-exl; |
| | 1199 | bins.fArray[i+1] = x+exh; |
| | 1200 | Double_t dxi = exh+exl; |
| | 1201 | if (i == 0) dx = dxi; |
| | 1202 | else if (dxi != dx) dx = 0; |
| | 1203 | } |
| | 1204 | TString name(Form("%s_%s", g1->GetName(), g2->GetName())); |
| | 1205 | TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle())); |
| | 1206 | TH1* h = 0; |
| | 1207 | if (dx != 0) { |
| | 1208 | h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]); |
| | 1209 | } |
| | 1210 | else { |
| | 1211 | h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray); |
| | 1212 | } |
| | 1213 | h->SetMarkerStyle(g2->GetMarkerStyle()); |
| | 1214 | h->SetMarkerColor(g1->GetMarkerColor()); |
| | 1215 | h->SetMarkerSize(0.9*g2->GetMarkerSize()); |
| | 1216 | |
| | 1217 | Double_t low = g2->GetX()[0]; |
| | 1218 | Double_t high = g2->GetX()[g2->GetN()-1]; |
| | 1219 | if (low > high) { Double_t t = low; low = high; high = t; } |
| | 1220 | for (Int_t i = 0; i < nBins; i++) { |
| | 1221 | Double_t x = g1->GetX()[i]; |
| | 1222 | if (x < low-dx || x > high+dx) continue; |
| | 1223 | Double_t c1 = g1->GetY()[i]; |
| | 1224 | Double_t e1 = g1->GetErrorY(i); |
| | 1225 | Double_t c2 = g2->Eval(x); |
| | 1226 | |
| | 1227 | h->SetBinContent(i+1, c1 / c2); |
| | 1228 | h->SetBinError(i+1, e1 / c2); |
| | 1229 | } |
| | 1230 | max = TMath::Max(RatioMax(h), max); |
| | 1231 | return h; |
| | 1232 | } |
| | 1233 | /* @} */ |
| | 1234 | //================================================================== |
| | 1235 | /** |
| | 1236 | * @{ |
| | 1237 | * @name Graphics utility functions |
| | 1238 | */ |
| | 1239 | /** |
| | 1240 | * Find an X axis in a pad |
| | 1241 | * |
| | 1242 | * @param p Pad |
| | 1243 | * @param name Histogram to find axis for |
| | 1244 | * |
| | 1245 | * @return Found axis or null |
| | 1246 | */ |
| | 1247 | TAxis* FindXAxis(TVirtualPad* p, const char* name) |
| | 1248 | { |
| | 1249 | TObject* o = p->GetListOfPrimitives()->FindObject(name); |
| | 1250 | if (!o) { |
| | 1251 | Warning("FindXAxis", "%s not found in pad", name); |
| | 1252 | return 0; |
| | 1253 | } |
| | 1254 | THStack* stack = dynamic_cast<THStack*>(o); |
| | 1255 | if (!stack) { |
| | 1256 | Warning("FindXAxis", "%s is not a THStack", name); |
| | 1257 | return 0; |
| | 1258 | } |
| | 1259 | if (!stack->GetHistogram()) { |
| | 1260 | Warning("FindXAxis", "%s has no histogram", name); |
| | 1261 | return 0; |
| | 1262 | } |
| | 1263 | TAxis* ret = stack->GetHistogram()->GetXaxis(); |
| | 1264 | return ret; |
| | 1265 | } |
| | 1266 | |
| | 1267 | //__________________________________________________________________ |
| | 1268 | /** |
| | 1269 | * Fix the apperance of the axis in a stack |
| | 1270 | * |
| | 1271 | * @param stack stack of histogram |
| | 1272 | * @param s Scaling factor |
| | 1273 | * @param ytitle Y axis title |
| | 1274 | * @param force Whether to draw the stack first or not |
| | 1275 | * @param ynDiv Divisions on Y axis |
| | 1276 | */ |
| | 1277 | void FixAxis(THStack* stack, Double_t s, const char* ytitle, |
| | 1278 | Int_t ynDiv=210, Bool_t force=true) |
| | 1279 | { |
| | 1280 | if (force) stack->Draw("nostack e1"); |
| | 1281 | |
| | 1282 | TH1* h = stack->GetHistogram(); |
| | 1283 | if (!h) return; |
| | 1284 | |
| | 1285 | h->SetXTitle("#eta"); |
| | 1286 | h->SetYTitle(ytitle); |
| | 1287 | TAxis* xa = h->GetXaxis(); |
| | 1288 | TAxis* ya = h->GetYaxis(); |
| | 1289 | if (xa) { |
| | 1290 | xa->SetTitle("#eta"); |
| | 1291 | // xa->SetTicks("+-"); |
| | 1292 | xa->SetTitleSize(s*xa->GetTitleSize()); |
| | 1293 | xa->SetLabelSize(s*xa->GetLabelSize()); |
| | 1294 | xa->SetTickLength(s*xa->GetTickLength()); |
| | 1295 | } |
| | 1296 | if (ya) { |
| | 1297 | ya->SetTitle(ytitle); |
| | 1298 | ya->SetDecimals(); |
| | 1299 | // ya->SetTicks("+-"); |
| | 1300 | ya->SetNdivisions(ynDiv); |
| | 1301 | ya->SetTitleSize(s*ya->GetTitleSize()); |
| | 1302 | ya->SetTitleOffset(ya->GetTitleOffset()/s); |
| | 1303 | ya->SetLabelSize(s*ya->GetLabelSize()); |
| | 1304 | } |
| | 1305 | } |
| | 1306 | /* @} */ |
| | 1307 | |
| | 1308 | |
| | 1309 | |
| | 1310 | //__________________________________________________________________ |
| | 1311 | Bool_t fShowOthers; // Show other data |
| | 1312 | Bool_t fShowAlice; // Show ALICE published data |
| | 1313 | Bool_t fShowRatios; // Show ratios |
| | 1314 | Bool_t fShowLeftRight;// Show asymmetry |
| | 1315 | UShort_t fRebin; // Rebinning factor |
| | 1316 | Bool_t fCutEdges; // Whether to cut edges |
| | 1317 | TString fTitle; // Title on plot |
| | 1318 | TString fHHDFile; // File name of old results |
| | 1319 | TNamed* fTrigString; // Trigger string (read, or set) |
| | 1320 | TNamed* fSNNString; // Energy string (read, or set) |
| | 1321 | TNamed* fSysString; // Collision system string (read or set) |
| | 1322 | TAxis* fVtxAxis; // Vertex cuts (read or set) |
| | 1323 | TH1* fForward; // Results |
| | 1324 | TH1* fForwardMC; // MC results |
| | 1325 | TH1* fForwardHHD; // Old results |
| | 1326 | TH1* fTruth; // MC truth |
| | 1327 | TH1* fCentral; // Central region data |
| | 1328 | TH1* fForwardSym; // Symmetric extension |
| | 1329 | TH1* fForwardMCSym; // Symmetric extension |
| | 1330 | TH1* fForwardHHDSym;// Symmetric extension |
| | 1331 | TH1* fTriggers; // Number of triggers |
| | 1332 | RangeParam* fRangeParam; // Parameter object for range zoom |
| | 1333 | |
| | 1334 | }; |
| | 1335 | |
| | 1336 | //=== Stuff for auto zooming ========================================= |
| | 1337 | void UpdateRange(dNdetaDrawer::RangeParam* p) |
| | 1338 | { |
| | 1339 | if (!p) { |
| | 1340 | Warning("UpdateRange", "No parameters %p", p); |
| | 1341 | return; |
| | 1342 | } |
| | 1343 | if (!p->fMasterAxis) { |
| | 1344 | Warning("UpdateRange", "No master axis %p", p->fMasterAxis); |
| | 1345 | return; |
| | 1346 | } |
| | 1347 | Int_t first = p->fMasterAxis->GetFirst(); |
| | 1348 | Int_t last = p->fMasterAxis->GetLast(); |
| | 1349 | Double_t x1 = p->fMasterAxis->GetBinCenter(first); |
| | 1350 | Double_t x2 = p->fMasterAxis->GetBinCenter(last); |
| | 1351 | //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2); |
| | 1352 | |
| | 1353 | if (p->fSlave1Axis) { |
| | 1354 | Int_t i1 = p->fSlave1Axis->FindBin(x1); |
| | 1355 | Int_t i2 = p->fSlave1Axis->FindBin(x2); |
| | 1356 | p->fSlave1Axis->SetRange(i1, i2); |
| | 1357 | p->fSlave1Pad->Modified(); |
| | 1358 | p->fSlave1Pad->Update(); |
| | 1359 | } |
| | 1360 | if (p->fSlave2Axis) { |
| | 1361 | Int_t i1 = p->fSlave2Axis->FindBin(x1); |
| | 1362 | Int_t i2 = p->fSlave2Axis->FindBin(x2); |
| | 1363 | p->fSlave2Axis->SetRange(i1, i2); |
| | 1364 | p->fSlave2Pad->Modified(); |
| | 1365 | p->fSlave2Pad->Update(); |
| | 1366 | } |
| | 1367 | TCanvas* c = gPad->GetCanvas(); |
| | 1368 | c->cd(); |
| | 1369 | } |
| | 1370 | |
| | 1371 | //____________________________________________________________________ |
| | 1372 | void RangeExec(dNdetaDrawer::RangeParam* p) |
| | 1373 | { |
| | 1374 | // Event types: |
| | 1375 | // 51: Mouse move |
| | 1376 | // 53: |
| | 1377 | // 1: Button down |
| | 1378 | // 21: Mouse drag |
| | 1379 | // 11: Mouse release |
| | 1380 | // dNdetaDrawer::RangeParam* p = |
| | 1381 | // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr); |
| | 1382 | Int_t event = gPad->GetEvent(); |
| | 1383 | TObject *select = gPad->GetSelected(); |
| | 1384 | if (event == 53) { |
| | 1385 | UpdateRange(p); |
| | 1386 | return; |
| | 1387 | } |
| | 1388 | if (event != 11 || !select || select != p->fMasterAxis) return; |
| | 1389 | UpdateRange(p); |
| | 1390 | } |
| | 1391 | |
| | 1392 | //=== Steering function ============================================== |
| | 1393 | void |
| | 1394 | DrawdNdeta(const char* filename="forward_dndeta.root", |
| | 1395 | Int_t flags=0xf, |
| | 1396 | const char* title="", |
| | 1397 | UShort_t rebin=5, |
| | 1398 | Bool_t cutEdges=false) |
| | 1399 | { |
| | 1400 | dNdetaDrawer d; |
| | 1401 | d.SetRebin(rebin); |
| | 1402 | d.SetCutEdges(cutEdges); |
| | 1403 | d.SetTitle(title); |
| | 1404 | d.SetHHDFile(""); |
| | 1405 | d.SetShowOthers(flags & 0x1); |
| | 1406 | d.SetShowAlice(flags & 0x2); |
| | 1407 | d.SetShowRatios(flags & 0x4); |
| | 1408 | d.SetShowLeftRight(flags & 0x8); |
| | 1409 | // Do the below if your input data does not contain these settings |
| | 1410 | // d.SetSNN(900); // Collision energy per nucleon pair (GeV) |
| | 1411 | // d.SetSys(1); // Collision system (1:pp, 2:PbPB) |
| | 1412 | // d.SetTrigger(1); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD) |
| | 1413 | // d.SetVertexRange(-10,10); // Collision vertex range (cm) |
| | 1414 | d.Run(filename); |
| | 1415 | } |
| | 1416 | //____________________________________________________________________ |
| | 1417 | // |
| | 1418 | // EOF |
| | 1419 | // |
| | 1420 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental