PWGLF/FORWARD/analysis2/scripts/SummaryMCCorrDrawer.C

   1 #include "SummaryDrawer.C"
   2
   3 /**
   4  * Class to draw a summary of the AOD production
   5  *
   6  * @par Input:
   7  * - The merged <tt>forward.root</tt> file.
   8  *   If the file isn't merged, it should still work.
   9  *
  10  * @par Output:
  11  * - A PDF file named after the input, but with <tt>.root</tt>
  12  *   replaced with <tt>pdf</tt>
  13  *
  14  */
  15 class SummaryMCCorrDrawer : public SummaryDrawer
  16 {
  17 public:
  18   enum EFlags {
  19     kEventInspector    = 0x001,
  20     kTrackDensity      = 0x002,
  21     kVertexBins        = 0x004,
  22     kResults           = 0x008,
  23     kCentral           = 0x010,
  24     kNormal            = 0x01F
  25   };
  26   SummaryMCCorrDrawer()
  27     : SummaryDrawer(),
  28       fSums(0),
  29       fResults(0)
  30   {}
  31   virtual ~SummaryMCCorrDrawer() {}
  32
  33   //__________________________________________________________________
  34   /**
  35    *
  36    *
  37    * @param fname
  38    * @param what
  39    */
  40   void Run(const char* fname, UShort_t what=kNormal)
  41   {
  42     // --- Open the file ---------------------------------------------
  43     TString filename(fname);
  44     TFile*  file = TFile::Open(filename, "READ");
  45     if (!file) {
  46       Error("Run", "Failed to open \"%s\"", filename.Data());
  47       return;
  48     }
  49
  50
  51     // --- Get top-level collection ----------------------------------
  52     fSums = GetCollection(file, "ForwardCorrSums");
  53     if (!fSums) return;
  54
  55     // --- Make our canvas -------------------------------------------
  56     TString pdfName(filename);
  57     pdfName.ReplaceAll(".root", ".pdf");
  58     CreateCanvas(pdfName, what & kLandscape);
  59
  60     // --- Make a Title page -------------------------------------------
  61     DrawTitlePage(file);
  62
  63     // --- Possibly make a chapter here ------------------------------
  64     if (what & kCentral && GetCollection(file, "CentralCorrSums"))
  65       MakeChapter("Forward");
  66
  67     // --- Set pause flag --------------------------------------------
  68     fPause = what & kPause;
  69
  70     // --- Do each sub-algorithm -------------------------------------
  71     if (what & kEventInspector) DrawEventInspector(fSums);
  72     if (what & kTrackDensity)   DrawTrackDensity(fSums);
  73     if (what & kVertexBins)     DrawVertexBins(true);
  74
  75     // --- Do the results ----------------------------------------------
  76     fResults = GetCollection(file, "ForwardCorrResults");
  77     if (!fResults) fResults = fSums; // Old-style
  78
  79     if (what & kResults) DrawResults(true);
  80
  81     // --- SPD clusters ----------------------------------------------
  82     if (what & kCentral) {
  83       // --- Get top-level collection --------------------------------
  84       fSums = GetCollection(file, "CentralCorrSums");
  85       if (fSums) {
  86         MakeChapter("Central");
  87         if (what & kEventInspector) DrawEventInspector(fSums);
  88         if (what & kTrackDensity)   DrawTrackDensity(fSums);
  89         if (what & kVertexBins)     DrawVertexBins(false);
  90       }
  91       else
  92         Warning("", "No CentralCorrSums found");
  93
  94       fResults = GetCollection(file, "CentralCorrResults");
  95       if (!fResults) fResults = fSums; // Old-style
  96       if (!fResults)
  97         Warning("", "No CentralCorrResults found");
  98
  99       if (what & kResults) DrawResults(false);
 100     }
 101
 102     CloseCanvas();
 103   }
 104 protected:
 105   //____________________________________________________________________
 106   void DrawTitlePage(TFile* file)
 107   {
 108     TCollection* c   = GetCollection(file, "ForwardCorrSums");
 109
 110     fBody->cd();
 111
 112     Double_t y = .9;
 113     TLatex* ltx = new TLatex(.5, y, "ESD+MC #rightarrow Corrections");
 114     ltx->SetTextSize(0.07);
 115     ltx->SetTextFont(62);
 116     ltx->SetTextAlign(22);
 117     ltx->SetNDC();
 118     ltx->Draw();
 119     y -= .075;
 120
 121     TCollection* ei = GetCollection(c, "fmdEventInspector");
 122     if (ei) {
 123       UShort_t sys=0;
 124       UShort_t sNN=0;
 125       Int_t field=0;
 126       ULong_t runNo=0;
 127       GetParameter(ei, "sys", sys);
 128       GetParameter(ei, "sNN", sNN);
 129       GetParameter(ei, "field", field);
 130       GetParameter(ei, "runNo", runNo);
 131
 132       TString tS; SysString(sys, tS); DrawParameter(y, "System", tS);
 133       TString tE; SNNString(sNN, tE); DrawParameter(y, "#sqrt{s_{NN}}", tE);
 134       DrawParameter(y, "L3 B field", Form("%+2dkG", field));
 135       DrawParameter(y, "Run #", Form("%6lu", runNo));
 136     }
 137
 138     PrintCanvas("MC Corrections");
 139   }
 140   //____________________________________________________________________
 141   TCollection* GetVertexList(TCollection* parent, const TAxis& axis, Int_t bin)
 142   {
 143     TString folder = TString::Format("vtx%+05.1f_%+05.1f",
 144                                      axis.GetBinLowEdge(bin),
 145                                      axis.GetBinUpEdge(bin));
 146     folder.ReplaceAll(".", "d");
 147     folder.ReplaceAll("+", "p");
 148     folder.ReplaceAll("-", "m");
 149
 150     TCollection* c = GetCollection(parent, folder);
 151     if (!c)
 152       Warning("DrawVertexBins", "List %s not found", folder.Data());
 153
 154     return c;
 155   }
 156   //____________________________________________________________________
 157   void DrawVertexBins(Bool_t forward)
 158   {
 159     // Info("DrawVertexBins", "Drawing %s vertex bins",
 160     //      forward ? "Forward" : "Central");
 161     TH1* vtxHist = GetH1(fSums, "vtxAxis");
 162     if (!vtxHist) return;
 163
 164     TAxis* vtxAxis = vtxHist->GetXaxis();
 165     for (Int_t i = 1; i <= vtxAxis->GetNbins(); i++) {
 166       // Info("", "Bin %d", i);
 167       TCollection* c = GetVertexList(fSums, *vtxAxis, i);
 168       if (!c) continue;
 169
 170       if (forward) {
 171         DivideForRings(true, true);
 172         for (UShort_t d = 1; d <= 3; d++) {
 173           for (UShort_t q = 0; q < (d == 1 ? 1 : 2); q++) {
 174             Char_t r = q == 0 ? 'I' : 'O';
 175             DrawInRingPad(d,r, GetH2(c, Form("FMD%d%c_cache",d,r)), "colz");
 176           }
 177         }
 178         DrawInPad(fBody, 2, GetH2(c, "primary"), "colz");
 179       }
 180       else {
 181         fBody->Divide(1,3, 0, 0);
 182
 183         DrawInPad(fBody, 1, GetH2(c, "hits"),     "colz");
 184         DrawInPad(fBody, 2, GetH2(c, "clusters"), "colz");
 185         DrawInPad(fBody, 3, GetH2(c, "primary"), "colz");
 186       }
 187       PrintCanvas(Form("%s sums - Vertex bin %+5.1f - %+5.1f",
 188                        (forward ? "Forward" : "Central"),
 189                        vtxAxis->GetBinLowEdge(i),
 190                        vtxAxis->GetBinUpEdge(i)));
 191     }
 192   }
 193   //____________________________________________________________________
 194   void DrawResults(Bool_t forward)
 195   {
 196     // Info("DrawVertexBins", "Drawing resulting %s vertex bins",
 197     //      forward ? "Forward" : "Central");
 198     TH1* vtxHist = GetH1(fSums, "vtxAxis");
 199     if (!vtxHist) return;
 200
 201     TAxis* vtxAxis = vtxHist->GetXaxis();
 202     for (Int_t i = 1; i <= vtxAxis->GetNbins(); i++) {
 203       TCollection* c = GetVertexList(fResults, *vtxAxis, i);
 204       if (!c) continue;
 205
 206       if (forward) {
 207         THStack* all = new THStack("all",
 208                                    "2^{nd} correction averaged over #phi");
 209         DivideForRings(true, true);
 210         for (UShort_t d = 1; d <= 3; d++) {
 211           for (UShort_t q = 0; q < (d == 1 ? 1 : 2); q++) {
 212             Char_t r = q == 0 ? 'I' : 'O';
 213             // TVirtualPad* p = RingPad(d, r);
 214             // p->cd();
 215             // p->Divide(1,2,0,0);
 216             TH2* h = GetH2(c, Form("FMD%d%c_vtxbin%03d",d,r,i));
 217             DrawInRingPad(d,r, h,"colz");
 218             // TVirtualPad* pp = p->cd(1);
 219             // TVirtualPad* ppp = p->cd(2);
 220             // ppp->SetRightMargin(pp->GetRightMargin());
 221             TH1* hh = h->ProjectionX();
 222             hh->Scale(1./ h->GetNbinsY());
 223             hh->SetFillColor(RingColor(d,r));
 224             hh->SetLineColor(RingColor(d,r));
 225             hh->SetMarkerColor(RingColor(d,r));
 226             hh->SetFillStyle(3001);
 227             hh->SetTitle(Form("#LT%s#GT", hh->GetTitle()));
 228             // DrawInPad(p,2, hh, "hist e");
 229             all->Add(hh, "hist e");
 230           }
 231           TVirtualPad* p = RingPad(0, '0');
 232           p->SetBottomMargin(0.10);
 233           p->SetLeftMargin(0.10);
 234           p->SetRightMargin(0.05);
 235           DrawInRingPad(0, 'O', all, "nostack");
 236         }
 237       }
 238       else {
 239         fBody->Divide(2,4,0,0);
 240         TH2* secMap               = GetH2(c, "secMap");
 241         TH2* secMapAlt            = GetH2(c, "secMapEff");
 242         if (!secMapAlt) secMapAlt = GetH2(c, "secMapHit");
 243         TH1* acc                  = GetH1(c, "acc");
 244         TH1* accAlt               = GetH1(c, "accEff");
 245         if (!accAlt) accAlt       = GetH1(c, "accHit");
 246
 247         TH1* secMapProj = secMap->ProjectionX();
 248         secMapProj->Scale(1./ secMap->GetNbinsY());
 249         secMapProj->Divide(acc);
 250         secMapProj->SetFillColor(kRed+1);
 251         secMapProj->SetFillStyle(3001);
 252         secMapProj->SetTitle(Form("#LT%s#GT/Acceptance", secMap->GetTitle()));
 253
 254         TH1* secMapAltProj = secMapAlt->ProjectionX();
 255         secMapAltProj->Scale(1./ secMapAlt->GetNbinsY());
 256         secMapAltProj->Divide(accAlt);
 257         secMapAltProj->SetFillColor(kBlue+1);
 258         secMapAltProj->SetFillStyle(3001);
 259         secMapAltProj->SetTitle(Form("#LT%s#GT/Acceptance",
 260                                      secMapAlt->GetTitle()));
 261
 262         Double_t secMapMax = TMath::Max(secMap->GetMaximum(),
 263                                         secMapAlt->GetMaximum());
 264         secMap->SetMaximum(secMapMax);
 265         secMapAlt->SetMaximum(secMapMax);
 266
 267         Double_t secMapProjMax = TMath::Max(secMapProj->GetMaximum(),
 268                                             secMapAltProj->GetMaximum());
 269         secMapProj->SetMaximum(secMapProjMax);
 270         secMapAltProj->SetMaximum(secMapProjMax);
 271
 272         acc->SetFillColor(kRed+1);
 273         acc->SetFillStyle(3001);
 274         accAlt->SetFillColor(kBlue+1);
 275         accAlt->SetFillStyle(3001);
 276
 277         Double_t accMax = TMath::Max(acc->GetMaximum(),accAlt->GetMaximum());
 278         acc->SetMaximum(accMax);
 279         accAlt->SetMaximum(accMax);
 280
 281         DrawInPad(fBody, 1, secMap,    "colz", kGridx);
 282         DrawInPad(fBody, 2, secMapAlt, "colz", kGridx);
 283
 284         TVirtualPad* p = fBody;
 285         TVirtualPad* pp = p->cd(1);
 286         TVirtualPad* ppp = p->cd(3);
 287         ppp->SetRightMargin(pp->GetRightMargin());
 288         DrawInPad(p,3, secMapProj, "hist", kGridx|kGridy);
 289
 290         ppp = p->cd(5);
 291         ppp->SetRightMargin(pp->GetRightMargin());
 292         DrawInPad(fBody, 5, acc, "", kGridx|kGridy);
 293
 294         pp = p->cd(2);
 295         pp->SetLeftMargin(0.10);
 296         ppp = p->cd(4);
 297         ppp->SetRightMargin(pp->GetRightMargin());
 298         ppp->SetLeftMargin(0.10);
 299         DrawInPad(p,4, secMapAltProj, "hist", kGridx|kGridy);
 300
 301         pp = p->cd(4);
 302         pp->SetLeftMargin(0.10);
 303         ppp = p->cd(6);
 304         ppp->SetRightMargin(pp->GetRightMargin());
 305         ppp->SetLeftMargin(0.10);
 306         DrawInPad(fBody, 6, accAlt, "", kGridx|kGridy);
 307
 308         DrawInPad(fBody, 7, GetH2(c, "diagnostics"), "colz");
 309
 310         TH2* ratio = static_cast<TH2*>(secMap->Clone("ratio"));
 311         ratio->Add(secMapAlt, -1);
 312         ratio->Divide(secMapAlt);
 313         ratio->SetTitle("Relative difference between maps");
 314         ratio->SetZTitle("#frac{S - S_{alt}}{S_{alt}}");
 315
 316         pp = p->cd(8);
 317         pp->SetLeftMargin(0.10);
 318         DrawInPad(fBody, 8, ratio, "colz");
 319       }
 320       PrintCanvas(Form("%s results - Vertex bin %+5.1f - %+5.1f",
 321                        (forward ? "Forward" : "Central"),
 322                        vtxAxis->GetBinLowEdge(i),
 323                        vtxAxis->GetBinUpEdge(i)));
 324     }
 325   }
 326   TCollection* fSums;
 327   TCollection* fResults;
 328 };
 329
 330 // #endif