PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C

   1 // ELossTimeTrain.C
   2 #ifndef NO_TRAIN
   3 #ifndef __CINT__
   4 #include <AliAnalysisManager.h>
   5 #include <fstream>
   6 #else
   7 class AliAnalysisManager;
   8 #endif
   9 #include "TrainSetup.C"
  10 #include "ParUtilities.C"
  11
  12 /**
  13  * Train to record time of each event
  14  *
  15  * @ingroup pwglf_forward_eventtime
  16  */
  17 class ELossTimeTrain : public TrainSetup
  18 {
  19 public:
  20   /**
  21    * Constructor
  22    *
  23    * @param name The name
  24    */
  25   ELossTimeTrain(const char* name="eventTime") : TrainSetup(name)
  26   {
  27     fOptions.Add("map", "FILE", "File containg map", "map.root");
  28     fOptions.Set("type", "ESD");
  29   }
  30   /**
  31    * Create our tasks
  32    */
  33   void CreateTasks(AliAnalysisManager*)
  34   {
  35     if (!fRailway->LoadLibrary("PWGLFforward2"))
  36       Fatal("CreateTasks", "Failed to load PWGLFforward2");
  37
  38     if (!ParUtilities::MakeScriptPAR(fRailway->Mode() == Railway::kLocal,
  39                                      "EventTimeTask.C",
  40                                      // Gui because of CDB - sigh!
  41                                      // XMLParser because of CDB
  42                                      "Gui,XMLParser,"
  43                                      "STEERBase,CDB,ESD,AOD,ANALYSIS,OADB,"
  44                                      "ANALYSISalice",
  45                                      fRailway))
  46       Fatal("","Failed to make support PAR");
  47     if (!fRailway->LoadLibrary("EventTimeTask"))
  48       Fatal("CreateTasks", "Failed to load EventTimeTask");
  49
  50     if (!ParUtilities::MakeScriptPAR(fRailway->Mode() == Railway::kLocal,
  51                                      "ELossTimeTask.C",
  52                                      "Gui,STEERBase,CDB,ESD,AOD,ANALYSIS,OADB,"
  53                                      "ANALYSISalice,PWGLFforward2,"
  54                                      "EventTimeTask",
  55                                      fRailway))
  56       Fatal("","Failed to make PAR");
  57     if (!fRailway->LoadLibrary("ELossTimeTask"))
  58       Fatal("CreateTasks", "Failed to load ELossTimeTask");
  59
  60     TString mapfile = fOptions.Get("map");
  61     gROOT->ProcessLine(Form("ELossTimeTask::Create(\"%s\")", mapfile.Data()));
  62
  63     fRailway->LoadAux(mapfile.Data(), true);
  64   }
  65   /**
  66    * Do not create a physics selection
  67    */
  68   // void CreatePhysicsSelection(Bool_t, AliAnalysisManager*) {}
  69   /**
  70    * Do not create a centrality selection
  71    */
  72   void CreateCentralitySelection(Bool_t) {}
  73   /**
  74    * Do not create an output handler
  75    */
  76   AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; }
  77   /**
  78    * The train class name
  79    *
  80    * @return Train class name
  81    */
  82   const char* ClassName() const { return "ELossTimeTrain"; }
  83   /**
  84    * Overloaded to create new dNdeta.C and dndeta.sh in the output
  85    * directory
  86    *
  87    * @param asShellScript
  88    */
  89   void SaveSetup(Bool_t asShellScript)
  90   {
  91     TrainSetup::SaveSetup(asShellScript);
  92
  93     SaveDraw();
  94   }
  95   void SaveDraw()
  96   {
  97     std::ofstream o("Draw.C");
  98     o << "// Written by " << ClassName() << "\n"
  99       << "void Draw(const char* fileName=\"AnalysisResults.root\")\n"
 100       << "{\n"
 101       << "  gSystem->AddIncludePath(\"-DNO_TRAIN -DSUMMARY\");\n"
 102       << "  const char* fwd = \"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
 103       << "  gSystem->AddIncludePath(Form(\"-I%s/scripts\", fwd));\n"
 104       << "  gROOT->SetMacroPath(Form(\"%s/trains:%s\", fwd,\n"
 105       << "                           gROOT->GetMacroPath()));\n"
 106       << "  gROOT->LoadMacro(\"ELossTimeTrain.C+\");\n"
 107       << "  ELossTimeSummary s;\n"
 108       << "  s.Run(fileName);\n"
 109       << "}\n"
 110       << std::endl;
 111     o.close();
 112   }
 113   void PostShellCode(std::ostream& f)
 114   {
 115     f << "  echo \"=== Draw results ...\"\n"
 116       << "  aliroot -l -b -q ${prefix}Draw.C\\(\\\"AnalysisResults.root\\\"\\)\n"
 117       << std::endl;
 118   }
 119 };
 120 #endif
 121 #ifdef SUMMARY
 122 # include <SummaryDrawer.C>
 123 # include <TColor.h>
 124
 125 /**
 126  * Draw summary of the above train
 127  *
 128  * @ingroup pwglf_forward_eventtime
 129  */
 130 struct ELossTimeSummary : public SummaryDrawer
 131 {
 132   enum EFlags {
 133     kEventInspector    = 0x001,
 134   };
 135   /**
 136    * Run the class
 137    *
 138    * @param fname Filename
 139    * @param flags Flags
 140    */
 141   void Run(const char* fname, UShort_t flags=0x01)
 142   {
 143     // --- Open the file -----------------------------------------------
 144     TString filename(fname);
 145     TFile* file = TFile::Open(filename.Data(), "READ");
 146     if (!file) {
 147       Error("Run", "Failed to open \"%s\"", filename.Data());
 148       return;
 149     }
 150     fPause         = flags & kPause;
 151
 152     // --- Make our canvas ---------------------------------------------
 153     TString pdfName(filename);
 154     pdfName.ReplaceAll(".root", ".pdf");
 155     CreateCanvas(pdfName, flags & kLandscape);
 156
 157     // --- Make title page -------------------------------------------
 158     TCollection* c = GetCollection(file, "elossTimeSums");
 159     DrawTitlePage(c);
 160
 161     if (flags & kEventInspector)    DrawEventInspector(c);
 162
 163     TH1* dt = GetH1(c, "dt");
 164     DrawInPad(fBody, 0, dt, "", kLogy);
 165     PrintCanvas("#Deltat");
 166
 167     const char* rings[] = { "FMD1i", "FMD2i", "FMD2o", "FMD3o", "FMD3i", 0 };
 168     const char** pring   = rings;
 169
 170     while (*pring) {
 171       DrawRing(c, *pring, dt);
 172       pring++;
 173     }
 174
 175     CloseCanvas();
 176   }
 177   /**
 178    * Draw the title page
 179    *
 180    * @param c Parent collection
 181    */
 182   void DrawTitlePage(TCollection* c)
 183   {
 184     fBody->cd();
 185
 186     Double_t y = .7;
 187     TLatex* ltx = new TLatex(.5, y, "#Deltat vs #Delta/#Delta_{mip}");
 188     ltx->SetTextSize(0.07);
 189     ltx->SetTextFont(62);
 190     ltx->SetTextAlign(22);
 191     ltx->SetNDC();
 192     ltx->Draw();
 193
 194     Double_t save = fParName->GetTextSize();
 195     fParName->SetTextSize(0.03);
 196     fParVal->SetTextSize(0.03);
 197     y = .6;
 198
 199     TCollection* fc = GetCollection(c, "fmdEventInspector");
 200     UShort_t sys;
 201     UShort_t sNN;
 202     ULong_t  runNo;
 203     GetParameter(fc, "sys",     sys);
 204     GetParameter(fc, "sNN",     sNN);
 205     GetParameter(fc, "runNo",   runNo);
 206
 207     DrawParameter(y, "Run #", Form("%lu", runNo));
 208     TString tS; SysString(sys, tS);      DrawParameter(y, "System", tS);
 209     TString tE; SNNString(sNN, tE);      DrawParameter(y, "#sqrt{s_{NN}}", tE);
 210
 211     PrintCanvas("Title page");
 212     fParName->SetTextSize(save);
 213     fParVal->SetTextSize(save);
 214   }
 215   /**
 216    * Draw a single ring
 217    *
 218    * @param c     Parent collection
 219    * @param ring  Ring name
 220    * @param dt    Histogram of delta time
 221    */
 222   void DrawRing(TCollection* c, const char* ring, TH1* dt)
 223   {
 224     TCollection* lring = GetCollection(c, ring);
 225     if (!lring) return;
 226
 227     TH2* h2 = GetH2(lring, "dtVsELoss");
 228     if (!h2) return;
 229
 230     THStack* stack = new THStack(ring, ring);
 231     // stack->SetTitle(ring);
 232
 233     THStack* ratios = new THStack(Form("Ratios for %s",ring), ring);
 234     // stack->SetTitle(ring);
 235
 236     Printf(ring);
 237     Int_t j = 2;
 238     TH1* first = 0;
 239     Double_t lfirst = 0;
 240     Double_t rmax = 0;
 241     Double_t max = 3;
 242     for (Int_t i = 1; i <= h2->GetNbinsY(); i++) {
 243       TH1*     h     = h2->ProjectionX(Form("%s_%03d", ring, i), i,i);
 244       Double_t logDt = h2->GetYaxis()->GetBinCenter(i);
 245
 246       Int_t    nFill = h->GetEntries();
 247       if (nFill <= 1000) continue;
 248       Double_t norm = dt->GetBinContent(i);
 249       if (norm <= 1e-6) {
 250         Warning("", "Normalization=%f<1e-6 but got "
 251                 "%d>1000 entries for log10(dt)=%5.3f", norm, nFill, logDt);
 252         continue;
 253       }
 254       if (!first && logDt > TMath::Log10(25.)) {
 255         lfirst = logDt;
 256         first = h;
 257       }
 258       // Info("", "Normalization is %f", norm);
 259       h->Sumw2();
 260       h->Scale(1. / norm);
 261       h->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt));
 262
 263       Float_t r, g, b;
 264       TColor::HSV2RGB((j-1)*45, 1, .8, r, g, b);
 265       Int_t col = TColor::GetColor(r, g, b);
 266       j++;
 267       h->SetLineColor(col);
 268       h->SetLineStyle(j % 3+1);
 269       h->SetLineWidth(2);
 270       // h->SetFillColor(col);
 271       // h->SetFillStyle(3002);
 272       stack->Add(h);
 273
 274       if (h == first) continue;
 275       TH1* rh = static_cast<TH1*>(h->Clone(Form("ratio%s", h->GetName())));
 276       // rh->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt));
 277       rh->Divide(first);
 278       for (Int_t k = 1; k <= rh->GetNbinsX(); k++) {
 279         if (rh->GetXaxis()->GetBinCenter(k) > max) break;
 280         rmax = TMath::Max(rmax, rh->GetBinContent(k));
 281       }
 282
 283       ratios->Add(rh);
 284     }
 285     Double_t savX = fParVal->GetX();
 286     Double_t savY = fParVal->GetY();
 287     fParVal->SetX(0.12);
 288     fParVal->SetY(0.12);
 289     fBody->Divide(1,2,0,0);
 290     DrawInPad(fBody,1,stack,"nostack hist", kLogy|kLegend);
 291     stack->GetXaxis()->SetRangeUser(-.1,max);
 292     stack->GetYaxis()->SetTitle("1/N_{ev} dN/d(#Delta/#Delta_{mip})");
 293
 294     fParVal->SetX(0.6);
 295     fParVal->SetY(0.4);
 296     DrawInPad(fBody,2,ratios,"nostack hist", kLegend);
 297     ratios->GetXaxis()->SetRangeUser(-.1, max);
 298     ratios->GetXaxis()->SetTitle("#Delta/#Delta_{mip}");
 299     ratios->GetYaxis()
 300       ->SetTitle(Form("X/(1/N_{ev}dN/d(#Delta/#Delta_{mip}))|_{%5.3f}",lfirst));
 301     Printf("Max: %f (%f)", ratios->GetMaximum(), rmax);
 302     ratios->SetMaximum(rmax*1.2);
 303
 304
 305     fParVal->SetX(savX);
 306     fParVal->SetY(savY);
 307     PrintCanvas(ring);
 308   }
 309 };
 310
 311 #endif
 312 // EOF