]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Added tasks and trains to investigate time-to-previous event.
authorChristian Holm Christensen <cholm@nbi.dk>
Mon, 9 Dec 2013 13:37:22 +0000 (14:37 +0100)
committerChristian Holm Christensen <cholm@nbi.dk>
Mon, 9 Dec 2013 13:37:22 +0000 (14:37 +0100)
Some fixes for the AOD and MultDist trains.

PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/trains/EventTimeTask.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C
PWGLF/FORWARD/analysis2/trains/MakeMultDistsTrain.C

diff --git a/PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C b/PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C
new file mode 100644 (file)
index 0000000..8a3f5b7
--- /dev/null
@@ -0,0 +1,372 @@
+#ifndef __CINT__
+# include <AliFMDEventInspector.h>
+# include <AliForwardCorrectionManager.h>
+# include <TH2.h>
+# include <AliESDFMD.h>
+# include <AliAODForwardMult.h>
+# include <AliESDEvent.h>
+# include <TFile.h>
+# include <TSystem.h>
+#else
+class AliFMDEventInspector;
+class TH2;
+class AliESDFMD;
+#endif
+#include <AliBaseESDTask.h>
+#include <AliForwardUtil.h>
+#define NO_TASK
+#define NO_SORTER
+#include "EventTimeTask.C"
+
+/**
+ * Task to analyse the energy loss in the FMD rings as function of the
+ * time to previous event.
+ * 
+ */
+struct ELossTimeTask : public AliBaseESDTask
+{
+  /** 
+   * Constructor - for I/O only
+   */
+  ELossTimeTask() 
+    : AliBaseESDTask(),
+      fEventInspector(),
+      fFMD1i(),
+      fFMD2i(),
+      fFMD2o(), 
+      fFMD3i(),
+      fFMD3o(),
+      fMap(0),
+      fDt(0)
+  {}
+  /** 
+   * Constructor 
+   * 
+   * @param name Name of task 
+   */     
+  ELossTimeTask(const char* name) 
+    : AliBaseESDTask(name, "ELossTimeTask", 
+                    &(AliForwardCorrectionManager::Instance())),
+      fEventInspector("fmdEventInspector"),
+      fFMD1i(1, 'i'),
+      fFMD2i(2, 'i'),
+      fFMD2o(2, 'o'), 
+      fFMD3i(3, 'i'),
+      fFMD3o(3, 'o'),
+      fMap(0),
+      fDt(0)
+  {
+  }
+  /** 
+   * Book output objects. Derived class should define this to book
+   * output objects on the processing output list @c fList before the
+   * actual event processing.  This is called on the master and on
+   * each slave.
+   * 
+   * If this member function returns false, the execution is stopped
+   * with a fatal signal.
+   *
+   * @return true on success. 
+   */
+  virtual Bool_t Book() 
+  {
+    fNeededCorrections = 0;
+    fExtraCorrections  = 0;
+
+    fDt = new TH1D("dt", "Time-to-last event (PS-triggered)", 60, 0, 15);
+    fDt->SetXTitle("log_{10}(#Deltat)");
+    fDt->SetFillColor(kYellow+2);
+    fDt->SetLineColor(kYellow+2);
+    fDt->SetFillStyle(3001);
+    fDt->SetDirectory(0);
+    fList->Add(fDt);
+
+    fFMD1i.Book(fList, fDt);
+    fFMD2i.Book(fList, fDt);
+    fFMD2o.Book(fList, fDt);
+    fFMD3i.Book(fList, fDt);
+    fFMD3o.Book(fList, fDt);
+
+    // Possibly re-read map
+    ReadMap("map.root");
+
+    return true;
+  }
+  /** 
+   * Process a single event
+   * 
+   * @param esd Input event 
+   * 
+   * @return true on success 
+   */
+  virtual Bool_t Event(AliESDEvent& esd)
+  {
+    Bool_t   lowFlux   = kFALSE;
+    UInt_t   triggers  = 0;
+    UShort_t ivz       = 0;
+    TVector3 ip;
+    Double_t cent      = 0;
+    UShort_t nClusters = 0;
+    UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux, 
+                                                ivz, ip, cent, nClusters);
+    if (found & AliFMDEventInspector::kNoEvent)    return false;
+    if (found & AliFMDEventInspector::kNoTriggers) return false;
+    if (found & AliFMDEventInspector::kNoSPD)      return false;
+    if (found & AliFMDEventInspector::kNoFMD)      return false;
+    if (found & AliFMDEventInspector::kNoVertex)   return false;
+    if (found & AliFMDEventInspector::kBadVertex)  return false;
+
+    // do not process pile-up, A, C, and E events 
+    if (triggers & AliAODForwardMult::kPileUp)     return false;
+    if (triggers & AliAODForwardMult::kA)          return false;
+    if (triggers & AliAODForwardMult::kC)          return false;
+    if (triggers & AliAODForwardMult::kE)          return false;
+  
+    // We want only the events found by off-line 
+    if (!(triggers & AliAODForwardMult::kOffline)) return false;
+    
+    // Perhaps we should also insist on MB only 
+    // if (fOnlyMB && (!(triggers & AliAODForwardMult::kInel))) return false;
+         
+    // Get FMD data 
+    AliESDFMD* esdFMD = esd.GetFMDData();
+    ULong64_t  period = esd.GetPeriodNumber();
+    ULong64_t  orbit  = esd.GetOrbitNumber();
+    ULong64_t  bc     = esd.GetBunchCrossNumber();
+    ULong64_t  full   = EventTimeData::EncodeFull(period, orbit, bc);
+    ULong64_t  dt     = fMap->Get(full);
+    Double_t   logDt  = TMath::Log10(25. * dt);
+    if (dt == EventTimeMap::kInvalidTime) {
+      logDt = 0;
+      Printf("!!! Failed to find dT for 0x%016llu", full);
+    }
+    // else 
+    //   Printf("=== 0x%016llu -> 0x%016llu", full, dt);
+    fDt->Fill(logDt);
+    
+    for (UShort_t d = 1; d <= 3; d++) { 
+      UShort_t nQ = d == 1 ? 1 : 2;
+      for (UShort_t q = 0; q < nQ; q++) { 
+       RingHistos* r = 0;
+       switch (d) { 
+       case 1: r = &fFMD1i; break;
+       case 2: r = (q == 0 ? &fFMD2i : &fFMD2o); break;
+       case 3: r = (q == 0 ? &fFMD3i : &fFMD3o); break;
+       }
+       r->Event(*esdFMD, logDt);
+      }
+    }
+    return true;
+  }
+  /** 
+   * Do the final analysis on the merged output. 
+   * 
+   * @return true on success
+   */
+  virtual Bool_t Finalize() 
+  { 
+    fDt = static_cast<TH1*>(fList->FindObject("dt"));
+
+    fFMD1i.Finalize(fList, fResults, fDt);
+    fFMD2i.Finalize(fList, fResults, fDt);
+    fFMD2o.Finalize(fList, fResults, fDt);
+    fFMD3i.Finalize(fList, fResults, fDt);
+    fFMD3o.Finalize(fList, fResults, fDt);
+    return true; 
+  }
+  /** 
+   * Get a reference to the event inspector. User must override this
+   * to return proper object
+   * 
+   * @return Reference to the event inspector 
+   */
+  virtual AliFMDEventInspector& GetEventInspector() { return fEventInspector; }
+  /** 
+   * Get a reference to the event inspector. User must override this
+   * to return proper object
+   * 
+   * @return Reference to the event inspector 
+   */
+  virtual const AliFMDEventInspector& GetEventInspector() const 
+  {
+    return fEventInspector;
+  }
+  /** 
+   * Read the map from timestamp to time-to-previous event 
+   * 
+   * @param filename File to read the map from 
+   * 
+   * @return true on success, false otherwise 
+   */
+  Bool_t ReadMap(const char* filename)
+  {
+    if (gSystem->AccessPathName(filename, kReadPermission)) {
+      // TSystem::AccessPathName returns false if file can be accessed!
+      Error("ReadMap", "File \"%s\" cannot be open for reading", filename);
+      return false;
+    }
+    Printf("Opening \"%s\" ...", filename);
+    TFile* file = TFile::Open(filename, "READ");
+    if (!file) { 
+      Error("ReadMap", "Failed to open file \"%s\"", filename);
+      return false;
+    }
+    Printf("Opened \"%s\" ...", filename);
+    TObject* o = file->Get("eventTimeMap");
+    if (!o) { 
+      Error("ReadMap", "Failed to get \"eventTimeMap\" from %s", filename);
+      return false;
+    }
+    Printf("Got object \"eventTimeMap\" from \"%s\" ...", filename);
+    if (!o->IsA()->InheritsFrom(EventTimeMap::Class())) { 
+      Error("ReadMap", "Object \"%s\" is not an EventTimeMap, but a %s", 
+           o->GetName(), o->ClassName());
+      return false;
+    }
+    Printf("Set the \"eventTimeMap\" to use");
+    if (fMap) { 
+      delete fMap;
+      fMap = 0;
+    }
+    fMap = static_cast<EventTimeMap*>(o);
+    file->Close();
+    return true;
+  }
+  /** 
+   * Create and connect the task 
+   * 
+   * @param mapfile File name of file containing timestamp map
+   * 
+   * @return true on connect
+   */
+  static Bool_t Create(const char* mapfile) 
+  {
+    ELossTimeTask* task = new ELossTimeTask("elossTime");
+    if (!task->ReadMap(mapfile)) return false;
+    task->Connect();
+    return true;
+  }
+protected:
+  /** 
+   * Dummy copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  ELossTimeTask(const ELossTimeTask& o) : AliBaseESDTask(o) {}
+  /** Our event inspector */
+  AliFMDEventInspector fEventInspector;
+  /** 
+   * Structure to hold per-ring histograms 
+   */
+  struct RingHistos : public AliForwardUtil::RingHistos
+  {
+    /** 
+     * Constructor - for I/O only
+     */
+    RingHistos() 
+      : AliForwardUtil::RingHistos(),
+       fDtVsELoss(0)
+    {}
+    /** 
+     * Constructor 
+     * 
+     * @param d detector number 
+     * @param r ring identifier 
+     */
+    RingHistos(UShort_t d, Char_t r) 
+      : AliForwardUtil::RingHistos(d,r),
+       fDtVsELoss(0)
+    {}
+    /** 
+     * Book histograms 
+     * 
+     * @param dir Parent list to add our list to
+     * 
+     * @return true on success 
+     */
+    Bool_t Book(TList* dir, TH1* dt)
+    {
+      TList* l = DefineOutputList(dir);
+      // Double_t dtBins[] = { 0, 1, 2e3, 5e5, 1.1e6, 5e6, 1e12, 1e20 };
+      fDtVsELoss = new TH2D("dtVsELoss", 
+                           Form("#Deltat vs #Delta/#Delta_{mip} - %s", 
+                                GetName()), 450, 0, 15, 
+                           dt->GetXaxis()->GetNbins(), 
+                           dt->GetXaxis()->GetXmin(), 
+                           dt->GetXaxis()->GetXmax());
+      fDtVsELoss->SetXTitle("#Delta/#Delta_{mip}");
+      fDtVsELoss->SetYTitle("log_{10}(#Deltat) [ns]");
+      fDtVsELoss->SetMarkerColor(Color());
+      fDtVsELoss->SetDirectory(0);
+      l->Add(fDtVsELoss);
+
+      return true;
+    }
+    /** 
+     * Process a single event
+     * 
+     * @param fmd FMD ESD data
+     * @param dt  Time-to-previous event
+     */
+    void Event(AliESDFMD& fmd, Double_t logDt)
+    {
+      const UShort_t nSec = NSector();
+      const UShort_t nStr = NStrip();
+      for (UShort_t s = 0; s < nSec; s++) { 
+       for (UShort_t t = 0; t < nStr; t++) { 
+         Float_t mult = fmd.Multiplicity(fDet,fRing,s,t);
+         if(mult == AliESDFMD::kInvalidMult || mult <= 0) 
+           continue;
+         fDtVsELoss->Fill(mult, logDt);
+       }
+      }
+    }
+    void Finalize(const TList* input, TList* output, TH1* dt)
+    {
+      TList* in  = GetOutputList(input);
+      TList* out = DefineOutputList(output);
+
+      fDtVsELoss = static_cast<TH2*>(in->FindObject("dtVsELoss"));
+
+      TH2* dtVsELoss = static_cast<TH2*>(fDtVsELoss->Clone());
+      dtVsELoss->SetDirectory(0);
+      dtVsELoss->SetZTitle("1/N_{events}");
+      dtVsELoss->Reset();
+      for (UShort_t j = 1; j <= dt->GetNbinsX(); j++) { 
+       Double_t norm = dt->GetBinContent(j);
+       if (norm <= 0) {
+         // Warning("Finalize", "Bin %d empty in dT", j);
+         continue;
+       }
+       for (UShort_t i = 1; i <= fDtVsELoss->GetNbinsX(); i++) { 
+         Double_t c = fDtVsELoss->GetBinContent(i, j);
+         Double_t e = fDtVsELoss->GetBinError(i, j);
+         dtVsELoss->SetBinContent(i,j,c/norm);
+         dtVsELoss->SetBinError(i,j,e/norm);
+       }
+      }
+      out->Add(dtVsELoss);
+    }
+    /** Our histogram */
+    TH2* fDtVsELoss;
+    ClassDef(RingHistos,1)
+  };
+  /** Container of FMD1i histograms */
+  RingHistos    fFMD1i;
+  /** Container of FMD2i histograms */
+  RingHistos    fFMD2i;
+  /** Container of FMD2o histograms */
+  RingHistos    fFMD2o;
+  /** Container of FMD3i histograms */
+  RingHistos    fFMD3i;
+  /** Container of FMD3o histograms */
+  RingHistos    fFMD3o;
+  /** Map from timestamp to time-to-previous event*/
+  EventTimeMap* fMap;
+  /** Distribution of log10(dt) */
+  TH1* fDt;
+  ClassDef(ELossTimeTask,2)
+};
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C b/PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C
new file mode 100644 (file)
index 0000000..2bfc467
--- /dev/null
@@ -0,0 +1,291 @@
+// ELossTimeTrain.C
+#ifndef NO_TRAIN
+#ifndef __CINT__
+#include <AliAnalysisManager.h>
+#include <fstream>
+#else
+class AliAnalysisManager;
+#endif
+#include "TrainSetup.C"
+#include "ParUtilities.C"
+
+/** 
+ * Train to record time of each event 
+ * 
+ */        
+class ELossTimeTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor 
+   * 
+   * @param name The name 
+   */
+  ELossTimeTrain(const char* name="eventTime") : TrainSetup(name)
+  { 
+    fOptions.Add("map", "FILE", "File containg map", "map.root");
+    fOptions.Set("type", "ESD");
+  }
+  /** 
+   * Create our tasks 
+   * 
+   * @param mgr Manager
+   */
+  void CreateTasks(AliAnalysisManager*)
+  {
+    if (!fHelper->LoadLibrary("PWGLFforward2")) 
+      Fatal("CreateTasks", "Failed to load PWGLFforward2");
+    
+    if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal,
+                                    "EventTimeTask.C",
+                                    // Gui because of CDB - sigh!
+                                    // XMLParser because of CDB 
+                                    "Gui,XMLParser,"
+                                    "STEERBase,CDB,ESD,AOD,ANALYSIS,OADB,"
+                                    "ANALYSISalice",
+                                    fHelper)) 
+      Fatal("","Failed to make support PAR");
+    if (!fHelper->LoadLibrary("EventTimeTask")) 
+      Fatal("CreateTasks", "Failed to load EventTimeTask");
+
+    if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal,
+                                    "ELossTimeTask.C",
+                                    "Gui,STEERBase,CDB,ESD,AOD,ANALYSIS,OADB,"
+                                    "ANALYSISalice,PWGLFforward2,"
+                                    "EventTimeTask",
+                                    fHelper)) 
+      Fatal("","Failed to make PAR");
+    if (!fHelper->LoadLibrary("ELossTimeTask")) 
+      Fatal("CreateTasks", "Failed to load ELossTimeTask");
+
+    TString mapfile = fOptions.Get("map");
+    gROOT->ProcessLine(Form("ELossTimeTask::Create(\"%s\")", mapfile.Data()));
+
+    fHelper->LoadAux(mapfile.Data(), true);
+  }
+  /** 
+   * Do not create a physics selection
+   */
+  // void CreatePhysicsSelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Do not create a centrality selection
+   */
+  void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Do not create an output handler
+   */
+  AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; }
+  /** 
+   * The train class name 
+   * 
+   * @return Train class name
+   */
+  const char* ClassName() const { return "ELossTimeTrain"; }
+  /** 
+   * Overloaded to create new dNdeta.C and dndeta.sh in the output 
+   * directory
+   * 
+   * @param asShellScript 
+   */
+  void SaveSetup(Bool_t asShellScript)
+  {
+    TrainSetup::SaveSetup(asShellScript);
+    
+    SaveDraw();
+  }
+  void SaveDraw()
+  {
+    std::ofstream o("Draw.C");
+    o << "// Written by " << ClassName() << "\n"
+      << "void Draw(const char* fileName=\"AnalysisResults.root\")\n"
+      << "{\n"
+      << "  gSystem->AddIncludePath(\"-DNO_TRAIN -DSUMMARY\");\n"
+      << "  const char* fwd = \"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
+      << "  gSystem->AddIncludePath(Form(\"-I%s/scripts\", fwd));\n"
+      << "  gROOT->SetMacroPath(Form(\"%s/trains:%s\", fwd,\n"
+      << "                           gROOT->GetMacroPath()));\n"
+      << "  gROOT->LoadMacro(\"ELossTimeTrain.C+\");\n"
+      << "  ELossTimeSummary s;\n"
+      << "  s.Run(fileName);\n"
+      << "}\n"
+      << std::endl;
+    o.close();
+  }
+  void PostShellCode(std::ostream& f)
+  {
+    f << "  echo \"=== Draw results ...\"\n"
+      << "  aliroot -l -b -q ${prefix}Draw.C\\(\\\"AnalysisResults.root\\\"\\)\n"
+      << std::endl;
+  }
+};
+#endif
+#ifdef SUMMARY
+# include <SummaryDrawer.C>
+# include <TColor.h>
+
+struct ELossTimeSummary : public SummaryDrawer 
+{
+  enum EFlags { 
+    kEventInspector    = 0x001, 
+  };
+
+  void Run(const char* fname, UShort_t flags=0x01)
+  {
+    // --- Open the file -----------------------------------------------
+    TString filename(fname);
+    TFile* file = TFile::Open(filename.Data(), "READ");
+    if (!file) { 
+      Error("Run", "Failed to open \"%s\"", filename.Data());
+      return;
+    }
+    fPause         = flags & kPause;
+    
+    // --- Make our canvas ---------------------------------------------
+    TString pdfName(filename);
+    pdfName.ReplaceAll(".root", ".pdf");
+    CreateCanvas(pdfName, flags & kLandscape);
+
+    // --- Make title page -------------------------------------------
+    TCollection* c = GetCollection(file, "elossTimeSums");
+    DrawTitlePage(c);
+
+    if (flags & kEventInspector)    DrawEventInspector(c);
+    
+    TH1* dt = GetH1(c, "dt");
+    DrawInPad(fBody, 0, dt, "", kLogy);
+    PrintCanvas("#Deltat");
+
+    const char* rings[] = { "FMD1i", "FMD2i", "FMD2o", "FMD3o", "FMD3i", 0 };
+    const char** pring   = rings;
+
+    while (*pring) {
+      DrawRing(c, *pring, dt);
+      pring++;
+    }
+
+    CloseCanvas();
+  }
+  void DrawTitlePage(TCollection* c)
+  {
+    fBody->cd();
+
+    Double_t y = .7;
+    TLatex* ltx = new TLatex(.5, y, "#Deltat vs #Delta/#Delta_{mip}");
+    ltx->SetTextSize(0.07);
+    ltx->SetTextFont(62);
+    ltx->SetTextAlign(22);
+    ltx->SetNDC();
+    ltx->Draw();
+
+    Double_t save = fParName->GetTextSize();
+    fParName->SetTextSize(0.03);
+    fParVal->SetTextSize(0.03);
+    y = .6;
+    
+    TCollection* fc = GetCollection(c, "fmdEventInspector");
+    UShort_t sys;
+    UShort_t sNN;
+    ULong_t  runNo;
+    GetParameter(fc, "sys",     sys);
+    GetParameter(fc, "sNN",     sNN);
+    GetParameter(fc, "runNo",   runNo);
+    
+    DrawParameter(y, "Run #", Form("%lu", runNo));
+    TString tS; SysString(sys, tS);      DrawParameter(y, "System", tS);
+    TString tE; SNNString(sNN, tE);      DrawParameter(y, "#sqrt{s_{NN}}", tE);
+
+    PrintCanvas("Title page");
+    fParName->SetTextSize(save);
+    fParVal->SetTextSize(save);
+  }
+  void DrawRing(TCollection* c, const char* ring, TH1* dt)
+  {
+    TCollection* lring = GetCollection(c, ring);
+    if (!lring) return;
+    
+    TH2* h2 = GetH2(lring, "dtVsELoss");
+    if (!h2) return;
+
+    THStack* stack = new THStack(ring, ring);
+    // stack->SetTitle(ring);
+
+    THStack* ratios = new THStack(Form("Ratios for %s",ring), ring);
+    // stack->SetTitle(ring);
+
+    Printf(ring);
+    Int_t j = 2;
+    TH1* first = 0;
+    Double_t lfirst = 0;
+    Double_t rmax = 0;
+    Double_t max = 3;
+    for (Int_t i = 1; i <= h2->GetNbinsY(); i++) { 
+      TH1*     h     = h2->ProjectionX(Form("%s_%03d", ring, i), i,i);
+      Double_t logDt = h2->GetYaxis()->GetBinCenter(i);
+
+      Int_t    nFill = h->GetEntries();
+      if (nFill <= 1000) continue;
+      Double_t norm = dt->GetBinContent(i);
+      if (norm <= 1e-6) {
+       Warning("", "Normalization=%f<1e-6 but got "
+               "%d>1000 entries for log10(dt)=%5.3f", norm, nFill, logDt);
+       continue;
+      }
+      if (!first && logDt > TMath::Log10(25.)) {
+       lfirst = logDt;
+       first = h;
+      }
+      // Info("", "Normalization is %f", norm);
+      h->Sumw2();
+      h->Scale(1. / norm);
+      h->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt));
+
+      Float_t r, g, b;
+      TColor::HSV2RGB((j-1)*45, 1, .8, r, g, b);
+      Int_t col = TColor::GetColor(r, g, b);
+      j++;
+      h->SetLineColor(col);
+      h->SetLineStyle(j % 3+1);
+      h->SetLineWidth(2);
+      // h->SetFillColor(col);
+      // h->SetFillStyle(3002);
+      stack->Add(h);
+
+      if (h == first) continue;
+      TH1* rh = static_cast<TH1*>(h->Clone(Form("ratio%s", h->GetName())));
+      // rh->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt));
+      rh->Divide(first);
+      for (Int_t k = 1; k <= rh->GetNbinsX(); k++) {
+       if (rh->GetXaxis()->GetBinCenter(k) > max) break;
+       rmax = TMath::Max(rmax, rh->GetBinContent(k));
+      }
+
+      ratios->Add(rh);
+    }
+    Double_t savX = fParVal->GetX();
+    Double_t savY = fParVal->GetY();
+    fParVal->SetX(0.12);
+    fParVal->SetY(0.12);
+    fBody->Divide(1,2,0,0);
+    DrawInPad(fBody,1,stack,"nostack hist", kLogy|kLegend);
+    stack->GetXaxis()->SetRangeUser(-.1,max);
+    stack->GetYaxis()->SetTitle("1/N_{ev} dN/d(#Delta/#Delta_{mip})");
+
+    fParVal->SetX(0.6);
+    fParVal->SetY(0.4);
+    DrawInPad(fBody,2,ratios,"nostack hist", kLegend);
+    ratios->GetXaxis()->SetRangeUser(-.1, max);
+    ratios->GetXaxis()->SetTitle("#Delta/#Delta_{mip}");
+    ratios->GetYaxis()
+      ->SetTitle(Form("X/(1/N_{ev}dN/d(#Delta/#Delta_{mip}))|_{%5.3f}",lfirst));
+    Printf("Max: %f (%f)", ratios->GetMaximum(), rmax);
+    ratios->SetMaximum(rmax*1.2);
+
+
+    fParVal->SetX(savX);
+    fParVal->SetY(savY);
+    PrintCanvas(ring);
+  }
+};
+
+#endif 
+// EOF
diff --git a/PWGLF/FORWARD/analysis2/trains/EventTimeTask.C b/PWGLF/FORWARD/analysis2/trains/EventTimeTask.C
new file mode 100644 (file)
index 0000000..95d49de
--- /dev/null
@@ -0,0 +1,776 @@
+// EventTimeTask.C
+#ifndef __CINT__
+# include <TTree.h>
+# include <TError.h>
+#else
+class TTree;
+// Force load of libGui
+class TGWindow;
+#endif
+class AliESDEvent;
+
+/**
+ * Data structure to be filled by task - one for each event
+ * 
+ */
+struct EventTimeData 
+{
+  /**
+   * Widths of field in the full timestamp 
+   */
+  enum { 
+    kBCWidth      = 12, 
+    kOrbitWidth   = 24, 
+    kPeriodWidth  = 28
+  };
+  /** Full time stamp */
+  ULong64_t fFull;
+  /** LDC creation time - standard time_t */
+  UInt_t    fTime;
+  /** Mask of detectors present in the event - FMD is 0x1000 */
+  UInt_t    fDetectors;
+  /** Type of event - 7 is physics */
+  UShort_t  fType;
+  /** 
+   * Create a branch in a tree 
+   * 
+   * @param tree Tree to create the branch in 
+   */
+  void CreateBranch(TTree* tree)
+  {
+    tree->Branch("event", &(this->fFull), 
+                "full/l:time/i:detector:type/s");
+  }
+  /** 
+   * Set the address of a branch for reading back objects from the tree
+   * 
+   * @param tree Tree to read back from 
+   */
+  void ReadBranch(TTree* tree) 
+  {
+    tree->SetBranchAddress("event", &(this->fFull));
+  }
+  /** 
+   * Utility function to encode the full time stamp from components. 
+   * 
+   * @param period Period counter (overflow of orbit counter)
+   * @param orbit  Orbit counter (period of 88us)
+   * @param bc     Bunch crossing number (period of 25ns)
+   * 
+   * @return Encoded full time stamp 
+   */
+  static ULong64_t EncodeFull(ULong64_t period, 
+                             ULong64_t orbit, 
+                             ULong64_t bc) 
+  {
+    const ULong64_t bcMask     = (1 << kBCWidth) - 1;
+    const ULong64_t orbitMask  = (1 << kOrbitWidth) - 1;
+    const ULong64_t periodMask = (1 << kPeriodWidth) - 1;
+    ULong64_t ret = ((bc     & bcMask)                  | 
+                    (orbit  & orbitMask)  <<  kBCWidth | 
+                    (period & periodMask) << (kBCWidth+kOrbitWidth));
+    return ret;
+  }
+  /** 
+   * Fill information from ESD into this data structure. 
+   * 
+   * @param esd  Event 
+   * @param dets List of active detectors in this event. 
+   */
+  void Fill(AliESDEvent* esd, UInt_t dets);
+};
+
+#ifndef NO_TASK
+# ifndef __CINT__
+#  include <AliESDEvent.h>
+# endif
+inline void EventTimeData::Fill(AliESDEvent* esd, UInt_t dets)
+{
+  ULong64_t period = esd->GetPeriodNumber();
+  ULong64_t orbit  = esd->GetOrbitNumber();
+  ULong64_t bc     = esd->GetBunchCrossNumber();
+  fType            = esd->GetEventType();
+  fTime            = esd->GetTimeStamp();//LDC time
+  fDetectors       = dets; // esd->GetDAQDetectorPattern();
+  fFull            = EncodeFull(period, orbit, bc);
+}
+# else
+inline void EventTimeData::Fill(AliESDEvent*, UInt_t) 
+{
+  Warning("Fill", "Calling empty method - shouldn't happen");
+}
+#endif
+
+#ifndef NO_TASK
+# ifndef __CINT__
+#  include <AliAnalysisManager.h>
+#  include <AliESDEvent.h>
+#  include <TTree.h>
+#  include <TH2.h>
+#  include <TList.h>
+#  include <AliCDBManager.h>
+#  include <AliCDBEntry.h>
+#  include <AliTriggerConfiguration.h>
+#  include <AliTriggerClass.h>
+#  include <AliTriggerCluster.h>
+#  include <AliDAQ.h>
+#  include <TObjArray.h>
+#  include <TDirectory.h>
+# else
+class AliAnalysisManager;
+class TTree;
+class AliTimeStamp;
+class TList;
+class TH2;
+# endif
+# include <AliAnalysisTaskSE.h>
+# include <vector>
+/**
+ * A task to record the unique timestamp of each event. 
+ *
+ * @par Input: ESD 
+ * @par Output:  A tree with a single branch 
+ */
+class EventTimeTask : public AliAnalysisTaskSE
+{
+public:
+  enum { 
+    kListSlot = 1,
+    kTreeSlot = 2
+  };
+  /** 
+   * Default CTOR - for I/O only.
+   */
+  EventTimeTask() 
+    : AliAnalysisTaskSE(), 
+      fTree(0),
+      fHistograms(0),
+      fDetVsType(0)
+  {
+  }
+  /** 
+   * User constructor 
+   * 
+   * @param name Name of task 
+   */
+  EventTimeTask(const char* name) 
+    : AliAnalysisTaskSE(name), 
+      fTree(0),
+      fHistograms(0),
+      fDetVsType(0)
+  {
+    DefineOutput(kListSlot, TList::Class());
+    DefineOutput(kTreeSlot, TTree::Class());
+    fBranchNames = "ESD:AliESDRun.,AliESDHeader.";
+  }
+  /** 
+   * Create user output objects.
+   *
+   * Called on each slave at start of job
+   */
+  void UserCreateOutputObjects()
+  {
+    Printf("Creating tree and histogram");
+    fHistograms = new TList();
+    fHistograms->SetOwner();
+    fHistograms->SetName("L");
+    
+    fDetVsType = new TH2D("detVsType", "Detector vs type", 
+                    16, -.5, 15.5, 31, -.5, 30.5);
+    fDetVsType->SetXTitle("Type");
+    fDetVsType->SetYTitle("Detector");
+    fDetVsType->SetDirectory(0);
+    fHistograms->Add(fDetVsType);
+    Printf("Histogram (%d,%f,%f)x(%d,%f,%f)",
+          fDetVsType->GetXaxis()->GetNbins(), 
+          fDetVsType->GetXaxis()->GetXmin(), 
+          fDetVsType->GetXaxis()->GetXmax(),
+          fDetVsType->GetYaxis()->GetNbins(), 
+          fDetVsType->GetYaxis()->GetXmin(), 
+          fDetVsType->GetYaxis()->GetXmax());
+
+    // TDirectory* savdir = gDirectory;
+    // Printf("Opening file at slot %d", kTreeSlot);
+    // OpenFile(kTreeSlot);
+    Printf("Make tree and disassociate from file");
+    fTree = new TTree("T", "T");
+    fTree->SetDirectory(0);
+    Printf("Create branch");
+    fData.CreateBranch(fTree);
+    // savdir->cd();
+    
+
+    PostData(kListSlot, fHistograms);
+    PostData(kTreeSlot, fTree);
+  }
+  /** 
+   * Analyse a single event 
+   * 
+   */
+  void UserExec(Option_t*)
+  {
+    static Bool_t first = true;
+    LoadBranches();
+
+    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+    if (!esd) return;
+    if (!esd->GetHeader()) return;     
+
+    if (first) {
+      LoadTriggerConfig(esd->GetRunNumber());
+      first = false;
+    }
+    ULong64_t mask = esd->GetTriggerMask();
+    UInt_t    dets = 0;
+    for (UShort_t i = 0; i < fDets.size(); i++) {
+      if ((1 << i) & mask) dets |= fDets[i];
+    }
+    // Printf("Event mask 0x%016llx -> 0x%08x", mask, dets);
+    fData.Fill(esd, dets);
+    fTree->Fill();
+
+    UInt_t type      = esd->GetEventType();
+    UInt_t detectors = dets;
+    for (UInt_t i = 0; i < 31; i++) { 
+      if ((1 << i) & detectors) fDetVsType->Fill(type, i);
+    }
+
+    PostData(kListSlot, fHistograms);
+    PostData(kTreeSlot, fTree);
+  }
+  void LoadTriggerConfig(Int_t runNo)
+  {
+    Printf("Loading trigger configuration for run %d",  runNo);
+    // --- Connect to CDB --------------------------------------------
+    AliCDBManager* cdb = AliCDBManager::Instance();
+    cdb->SetDefaultStorageFromRun(runNo);
+    cdb->SetRun(runNo);
+    
+    // --- Get entry -------------------------------------------------
+    AliCDBEntry* entry = cdb->Get("GRP/CTP/Config");
+    if (!entry || !entry->GetObject()) { 
+      Warning("LoadTriggerConfig", "Couldn't get trigger configuration");
+      return;
+    }
+    AliTriggerConfiguration* config = 
+      static_cast<AliTriggerConfiguration*>(entry->GetObject());
+
+    // --- Get the classes, and resize cache -------------------------
+    const TObjArray& clss = config->GetClasses();
+    fDets.resize(clss.GetEntries());
+
+    // --- Loop over configurations ----------------------------------
+    TIter next(&clss);
+    AliTriggerClass* cls = 0;
+    while ((cls = static_cast<AliTriggerClass*>(next()))) {
+      Int_t              mask    = cls->GetMask();
+      AliTriggerCluster* cluster = cls->GetCluster();
+      if (!cluster) { 
+       Warning("LoadTriggerConfig", 
+               "Failed to get trigger cluster for %s", cls->GetName());
+       continue;
+      }
+      TString names = cluster->GetDetectorsInCluster();
+      if (names.IsNull()) { 
+       Warning("LoadTriggerConfig", "No detectors for cluster %s",
+               cls->GetName());
+       continue;
+      }
+      UInt_t dets = AliDAQ::DetectorPattern(names);
+      UInt_t id   = UInt_t(TMath::Log2(mask));
+      fDets[id]   = dets;
+      Printf("Trigger mask 0x%08x (%3d): 0x%08x (%s)", 
+            mask, id, dets, names.Data());
+    }
+    // cdb->Destroy();
+  }
+  /** 
+   * Register with manager and connect output containers
+   * 
+   * @param mgr Manager 
+   */
+  void Connect(AliAnalysisManager* mgr)
+  {
+    mgr->AddTask(this);
+    mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
+    AliAnalysisDataContainer* contTree = 
+      mgr->CreateContainer("T", TTree::Class(), 
+                          AliAnalysisManager::kOutputContainer,
+                          "time.root");
+    AliAnalysisDataContainer* contList = 
+      mgr->CreateContainer("L", TList::Class(), 
+                          AliAnalysisManager::kOutputContainer,
+                          "hist.root");
+    mgr->ConnectOutput(this, kTreeSlot, contTree);
+    mgr->ConnectOutput(this, kListSlot, contList);
+  }
+  /** 
+   * Create an instance of this task, and register it and it's
+   * outputs.
+   */
+  static void Create()
+  {
+    EventTimeTask* task = new EventTimeTask("eventTime");
+    task->Connect(AliAnalysisManager::GetAnalysisManager());
+  }
+  TTree*              fTree;       // Our tree
+  EventTimeData       fData;       // Our data 
+  TList*              fHistograms; // List
+  TH2D*               fDetVsType;  // Histogram
+  std::vector<UInt_t> fDets;       // Per-trigger-bit detector mask 
+  
+  ClassDef(EventTimeTask,3);
+};
+#endif // NO_TASK 
+
+#ifndef NO_MAP
+#include <utility>
+#include <map>
+
+typedef std::pair<ULong64_t,ULong64_t> EventTimeMapPair;
+
+/**
+ * A map of event time-stamp to distance to previous event
+ * 
+ */
+struct EventTimeMap : public TObject
+{
+  /** Map type */
+  typedef std::map<ULong64_t,ULong64_t> Map;
+  /** Key type */
+  typedef Map::key_type Key;
+  /** Mapped value type */
+  typedef Map::mapped_type Value;
+  /** Element type */
+  typedef Map::value_type Pair;
+  /** Iterator type */
+  typedef Map::iterator   Iterator;
+  /** Constant iterator type */
+  typedef Map::const_iterator ConstIterator;
+  /** 
+   * Constructor
+   */
+  EventTimeMap() : TObject(), fMap() {}
+  /** 
+   * Destructor 
+   */
+  virtual ~EventTimeMap() {}
+  /** 
+   * Get name of this object - always the same 
+   * 
+   * @return The string "eventTimeMap"
+   */
+  const char* GetName() const { return "eventTimeMap"; }
+  /** 
+   * Element access.  If the key @a k doesn't already exist, it is
+   * created
+   * 
+   * @param k Key 
+   * 
+   * @return A pair of key and value 
+   */
+  Value& operator[](const Key& k)
+  {
+    return fMap[k];
+  } 
+  /** 
+   * Find the element whos key is @a k 
+   * 
+   * @param k Key to look for
+   * 
+   * @return Iterator pointing to element, or end of container if not found 
+   */
+  Iterator Find(const Key& k)
+  {
+    return fMap.find(k);
+  }
+  /** 
+   * Find the element whos key is @a k 
+   * 
+   * @param k Key to look for
+   * 
+   * @return Iterator pointing to element, or end of container if not found 
+   */
+  ConstIterator Find(const Key& k) const
+  {
+    return fMap.find(k);
+  }
+  /** 
+   * Get forward iterator pointing beginning of the container
+   *
+   * @return Iterator to start of container
+   */
+  Iterator Begin()
+  {
+    return fMap.begin();
+  }
+  /** 
+   * Get forward iterator pointing beginning of the container
+   *
+   * @return Iterator to start of container
+   */
+  ConstIterator Begin() const
+  {
+    return fMap.begin();
+  }
+  /** 
+   * Get forward iterator pointing just beyond the end of the container
+   *
+   * @return Iterator just beyond container
+   */
+  Iterator End()
+  {
+    return fMap.end();
+  }
+  /** 
+   * Get forward iterator pointing just beyond the end of the container
+   *
+   * @return Iterator just beyond container
+   */
+  ConstIterator End() const
+  {
+    return fMap.end();
+  }
+  enum { 
+    kInvalidTime = 0xFFFFFFFFFFFFFFFF
+  };
+  /**
+   * Get the time difference to previous event from a event with a
+   * given time stamp.
+   *
+   * @param timestamp Time stamp of the event 
+   *
+   * @return time difference or kInvalidTime if not found 
+   */
+  Value Get(const Key& timestamp) const 
+  {
+    ConstIterator i = Find(timestamp);
+    if (i == End()) return kInvalidTime;
+    return i->second;
+  }
+  /** 
+   * Set the time to previous event for a given event time stamp
+   * 
+   * @param timestamp Event time stamp 
+   * @param diff      Time to previous event 
+   */
+  void Set(const Key& timestamp, const Value& diff)
+  {
+    this->operator[](timestamp) = diff;
+  }
+  ULong64_t Size() const { return fMap.size(); }
+  /** Our embeded map */
+  Map fMap;
+  ClassDef(EventTimeMap,1)
+};
+#endif // NO_MAP
+
+#ifndef NO_SORTER
+# ifndef __CINT__
+#  include <TFile.h>
+#  include <TArrayL64.h>
+#  include <TMath.h>
+#  include <TParameter.h>
+#  include <TCanvas.h>
+#  include <TH1.h>
+# else 
+class TFile;
+class TTree;
+class TCanvas;
+# endif
+
+/** 
+ * A class to sort the tree and generate our timestamp->dT map.
+ */
+struct EventTimeSorter 
+{
+  TTree*        fTree; 
+  EventTimeData fData;
+  ULong64_t     fMin;
+  ULong64_t     fMax;
+
+  /** 
+   * Constructor 
+   */
+  EventTimeSorter() : fTree(0), fData(), fMin(0xFFFFFFFFFFFFFFFF), fMax(0) {}
+  /** 
+   * Progress meter 
+   * 
+   * @param cur    Current entry
+   * @param total  Total number of entries
+   */
+  void Progress(Long64_t cur, Long64_t total) const
+  {
+    static UShort_t old = 0;
+    if (cur == 0) old = 0;
+    UShort_t now = UShort_t(100 * (cur + 1 == total ? 1 : 
+                                  Double_t(cur)/total));
+    if (now > old) 
+      std::cout << "\rLooping over " << total << " entries: ... "
+               << now << '%' << std::flush;
+    if (now >= 100) std::cout << std::endl;
+    old = now;
+  }
+  /** 
+   * Connect a tree 
+   * 
+   * @param inputName 
+   * @param treeName 
+   * 
+   * @return 
+   */
+  Bool_t OpenInput(const char* inputName, const char* treeName)
+  {
+    CloseInput();
+
+    // --- Get input -------------------------------------------------
+    TFile* input = TFile::Open(inputName, "READ");
+    if (!input) { 
+      Error("Run", "Failed to open \"%s\"", inputName);
+      return false;
+    }
+    
+    fTree = static_cast<TTree*>(input->Get(treeName));
+    if (!fTree) { 
+      Error("Run", "Couldn't get tree \"%s\" from \"%s\"", treeName,inputName);
+      return false;
+    }
+
+    // --- Set branch address ---------------------------------------
+    fData.ReadBranch(fTree);
+
+    return true;
+  }
+  /** 
+   * Disconnect tree
+   */          
+  void CloseInput()
+  {
+    if (!fTree) return;
+    TFile* file = fTree->GetCurrentFile();
+    if (file) file->Close();
+    fTree = 0;
+  }
+  /** 
+   * Run the sorter 
+   * 
+   * @param inputName   Input file name 
+   * @param outputName  Output file name 
+   * @param treeName    Tree name 
+   * 
+   * @return true on success
+   */
+  Bool_t Run(const char* inputName, const char* outputName, 
+            const char* treeName="T") 
+  {
+    if (!OpenInput(inputName, treeName)) return false;
+
+    // --- Loop over the data ----------------------------------------
+    Long64_t   nEnt    = fTree->GetEntries();
+    Long64_t   n       = 0;
+    ULong64_t* values  = new ULong64_t[nEnt];
+    UInt_t*    dets    = new UInt_t[nEnt];
+    UShort_t*  types   = new UShort_t[nEnt];
+    UInt_t*    times   = new UInt_t[nEnt];
+    for (Long64_t iEnt = 0; iEnt < nEnt; iEnt++) { 
+      Progress(iEnt, nEnt);
+
+      fTree->GetEntry(iEnt);
+
+      if (!(fData.fDetectors & 0x1000)) 
+       // No FMD read-out 
+       continue;
+      if (fData.fType != 7) { 
+       // Ignore all but physics events 
+       Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType);
+       // continue;
+      }
+
+      // Store values 
+      Int_t j = n;
+      values[j] = fData.fFull;
+      dets[j]   = fData.fDetectors;
+      types[j]  = fData.fType;
+      times[j]  = fData.fTime;
+      n++;
+    }
+    
+    // --- Now sort all values ---------------------------------------
+    TArrayL64 index(n);
+    TMath::Sort(n, values, index.fArray, false);
+    
+    // Maps the unique time to the distance to previous event 
+    EventTimeMap  eventTimeMap;
+    ULong64_t     key   = values[index[0]];
+    ULong64_t     prev  = 0;
+    ULong64_t     dt    = key-prev;
+    ULong64_t     min   = 0xFFFFFFFFFFFFFFFF;
+    ULong64_t     max   = 0;
+    eventTimeMap[key]   = dt; 
+    for (Long64_t i = 1; i < n; i++) { 
+      Long64_t j        = index[i];
+      Long64_t k        = index[i-1];
+      key               = values[j];
+      prev              = values[k];
+      dt                = key - prev;
+      if (dt <= 0) {
+       Printf("0x%016llx==0x%016llx -> dt=%10llu [%10lld %10lld] "
+              "(det: 0x%08x 0x%08x  type: 0x%04x 0x%04x  time: %08d %08d)",
+              key, prev, dt, j, k, dets[j], dets[k], 
+              types[j], types[k], times[j], times[k]);
+       // continue;
+      }
+      eventTimeMap[key] = dt;
+      min               = TMath::Min(dt, min);
+      max               = TMath::Max(dt, max);
+    }
+    std::cout << "Range is [" << min << "," << max << ']' << std::endl;
+
+    TFile* output = TFile::Open(outputName, "RECREATE");
+    if (!output) { 
+      Error("Run", "Failed to create output file \"%s\"", outputName);
+      return false;
+    }
+    fMin = min;
+    fMax = max;
+    eventTimeMap.Write();
+    output->Write();
+    output->Close();
+    
+    delete [] values; 
+
+    CloseInput();
+
+    return true;
+  }
+  Bool_t Test(const char* inputName, const char* outputName, 
+             const char* treeName="T") 
+  {
+    if (!OpenInput(inputName, treeName)) return false;
+
+    // --- Get our map -----------------------------------------------
+    TFile* output = TFile::Open(outputName, "UPDATE");
+    if (!output) { 
+      Error("Test", "Failed to open \"%s\"", outputName);
+      return false;
+    }
+
+    EventTimeMap* eventTimeMap = 
+      static_cast<EventTimeMap*>(output->Get("eventTimeMap"));
+    if (!eventTimeMap) { 
+      Error("Test", "Failed to get \"%s\" from \"%s\"",
+           "eventTimeMap", outputName);
+      return false;
+    }
+
+    // --- Histograms --------------------------------------------------
+    ULong64_t mmin = TMath::Min(25*fMin, 900000ULL);
+    ULong64_t mmax = TMath::Min(25*fMax, 110000000ULL);
+    TH1*      diff = new TH1D("diff", "Time-to-last-event (10#mus bins w/FMD)", 
+                             (mmax-mmin)/10000, mmin, mmax);
+    diff->SetStats(0);
+    diff->SetXTitle("#Deltat [ns]");
+    diff->SetFillColor(kRed+1);
+    diff->SetFillStyle(3001);
+    diff->SetDirectory(0);
+    
+
+    TH1* ldiff = new TH1D("ldiff", "log(#Deltat) - Events w/FMD",
+                         300, 0, 15);
+    ldiff->SetStats(0);
+    ldiff->SetXTitle("log_{10}(#Deltat) [ns]");
+    ldiff->SetFillColor(kGreen+1);
+    ldiff->SetFillStyle(3001);
+    ldiff->SetDirectory(0);
+                        
+    // --- Loop over the data ----------------------------------------
+    Long64_t   nEnt    = fTree->GetEntries();
+    Long64_t   nZero   = 0;
+    Long64_t   nMiss   = 0;
+    Long64_t   nGood   = 0;
+    Long64_t   nNoFMD  = 0;
+    for (Long64_t iEnt = 0; iEnt < /*10*/ nEnt; iEnt++) { 
+      Progress(iEnt, nEnt);
+      fTree->GetEntry(iEnt);
+      
+      if (!(fData.fDetectors & 0x1000)) {
+       // No FMD read-out 
+       nNoFMD++;
+       continue;
+      }
+      if (fData.fType != 7) { 
+       // Ignore all but physics events 
+       Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType);
+       // continue;
+      }
+      
+      // Look-up the timestamp 
+      ULong64_t value = fData.fFull;
+      ULong64_t dT    = eventTimeMap->Get(value);
+      if (dT == EventTimeMap::kInvalidTime) {
+       Warning("Test", "Value %llu not found", value);
+       ldiff->Fill(1);
+       nMiss++;
+       continue;
+      }
+      if (dT <= 0) { 
+#if 0
+       Warning("Test", "Impossible dt=%llu found for %llu (0x%0x %2d)", 
+               dT, value, fData.fDetectors, fData.fType);
+#endif
+       ldiff->Fill(1);
+       nZero++;
+       continue;
+      }
+      nGood++;
+      Double_t  dt    = 25.*dT;
+      diff->Fill(dt);
+      Double_t  logDt = TMath::Log10(dt);
+      ldiff->Fill(logDt);
+    }
+    CloseInput();
+    Printf("missing: %llu, bad: %llu, good: %llu, no FMD: %llu, all: %llu", 
+          nMiss,nZero,nGood,nNoFMD,nEnt);
+
+    // --- Draw results ----------------------------------------------
+    TCanvas* c = new TCanvas("c", "C");
+    c->SetTopMargin(0.01);
+    c->SetRightMargin(0.01);
+    c->Divide(2,1); // ,0,0);
+    TVirtualPad* p = c->cd(2);
+    // p->SetRightMargin(0.10);
+    p->SetLogy();
+    ldiff->DrawCopy();
+    
+    p = c->cd(1);
+    p->SetLogy();
+    p->SetLogx();
+    diff->DrawCopy();
+
+    c->Print("dt.png");
+    c->Print("dt.C");
+
+    // --- Write histogram to file -----------------------------------
+    output->cd();
+    diff->Write();
+    output->Write();
+    output->Close();
+    
+    return true;
+  }
+};
+#endif // NO_SORTER 
+#ifdef __MAKECINT__
+# ifndef NO_MAP
+#  pragma link C++ class std::pair<ULong64_t,ULong64_t>+;
+#  pragma link C++ class std::map<ULong64_t,ULong64_t>+;
+# endif
+# ifndef NO_TASK
+#  pragma link C++ class std::vector<UInt_t>+;
+# endif
+#endif
+
+// EOF
+        
+        
diff --git a/PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C b/PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C
new file mode 100644 (file)
index 0000000..1c9f75a
--- /dev/null
@@ -0,0 +1,104 @@
+// EventTimeTrain.C
+#ifndef __CINT__
+#include <AliAnalysisManager.h>
+#include <fstream>
+#else
+class AliAnalysisManager;
+#endif
+#include "TrainSetup.C"
+#include "ParUtilities.C"
+
+/** 
+ * Train to record time of each event 
+ * 
+ */        
+class EventTimeTrain : public TrainSetup
+{
+public:
+  /** 
+   * Constructor 
+   * 
+   * @param name The name 
+   */
+  EventTimeTrain(const char* name="eventTime") : TrainSetup(name)
+  { 
+    fOptions.Set("type", "ESD");
+  }
+  /** 
+   * Create our tasks 
+   * 
+   * @param mgr Manager
+   */
+  void CreateTasks(AliAnalysisManager*)
+  {
+    if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal,
+                                    "EventTimeTask.C",
+                                    // Gui because of CDB 
+                                    // XMLParser because of CDB
+                                    // CDB because look-up of trigger config
+                                    "Gui,XMLParser,"
+                                    "STEERBase,CDB,ESD,AOD,ANALYSIS,OADB,"
+                                    "ANALYSISalice",
+                                    fHelper)) 
+      Fatal("","Failed to make PAR");
+    fHelper->LoadLibrary("EventTimeTask");
+    gROOT->ProcessLine("EventTimeTask::Create()");
+  }
+  /** 
+   * Do not create a physics selection
+   */
+  void CreatePhysicsSelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Do not create a centrality selection
+   */
+  void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {}
+  /** 
+   * Do not create an output handler
+   */
+  AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; }
+  /** 
+   * The train class name 
+   * 
+   * @return Train class name
+   */
+  const char* ClassName() const { return "EventTimeTrain"; }
+  /** 
+   * Overloaded to create new dNdeta.C and dndeta.sh in the output 
+   * directory
+   * 
+   * @param asShellScript 
+   */
+  void SaveSetup(Bool_t asShellScript)
+  {
+    TrainSetup::SaveSetup(asShellScript);
+    
+    SaveSort();
+  }
+  void SaveSort()
+  {
+    std::ofstream o("Sort.C");
+    o << "// Written by " << ClassName() << "\n"
+      << "void Sort(const char* prefix=\"\",\n"
+      << "          const char* fileName=\"time.root\",\n"
+      << "          const char* outName=\"map.root\",\n"
+      << "          const char* treeName=\"T\")\n"
+      << "{\n"
+      << "  gSystem->AddIncludePath(\"-DNO_TASK -I$ALICE_ROOT/include\");\n"
+      << "  TString mac(\"EventTimeTask/EventTimeTask.C+g\");\n"
+      << "  if (prefix && prefix[0] != '\\0') mac.Prepend(prefix);\n"
+      << "  gROOT->LoadMacro(mac);\n"
+      << "  EventTimeSorter s;\n"
+      << "  if (!s.Run(fileName,outName,treeName)) return;\n"
+      << "  s.Test(fileName,outName,treeName);\n"
+      << "}\n"
+      << std::endl;
+    o.close();
+  }
+  void PostShellCode(std::ostream& f)
+  {
+    f << "  echo \"=== Sort results ...\"\n"
+      << "  aliroot -l -b -q ${prefix}Sort.C\\(\\\"${prefix}\\\"\\)\n"
+      << std::endl;
+  }
+};
+// EOF
index 3f97b380e4fa69faff8fa1d0befa275617097027..56ee98b7c3c170ead137a5bbca6c507895eeff66 100644 (file)
@@ -264,7 +264,7 @@ protected:
       << "//   0x200     Pause\n"
       << "//\n"
       << "void Summarize(const char* filename=\"forward.root\",\n"
-      << "               UShort_t what=0x0FF)\n"
+      << "               UShort_t what=0x1FF)\n"
       << "{\n"
       << "  const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
       << "  gROOT->LoadMacro(Form(\"%s/DrawAODSummary.C\",fwd));\n"
index f8d702b427e4c58e00e98ce014db40493a83330e..60532c56948c88b6525122f2ed94200318490521 100644 (file)
@@ -33,7 +33,7 @@ public:
     fOptions.Add("trig",       "TYPE", "Trigger type", "V0AND");
     fOptions.Add("vzMin",      "CENTIMETER", "Min Ip Z", -4);
     fOptions.Add("vzMax",      "CENTIMETER", "Max Ip Z", +4);
-    fOptions.Add("phi-acc",    "Use stored phi acceptance", true);
+    fOptions.Add("phi-acc",    "Use stored phi acceptance", false);
     fOptions.Add("asymmetric", "Make asymmetric (+/-) bins", false);
   }
 protected: