PWGLF/FORWARD/analysis2/trains/EventTimeTask.C

   1 // EventTimeTask.C
   2 #ifndef __CINT__
   3 # include <TTree.h>
   4 # include <TError.h>
   5 # include <TChain.h>
   6 #else
   7 class TTree;
   8 // Force load of libGui
   9 class TGWindow;
  10 #endif
  11 class AliESDEvent;
  12
  13 /**
  14  * Data structure to be filled by task - one for each event
  15  *
  16  * @ingroup pwglf_forward_eventtime
  17  */
  18 struct EventTimeData
  19 {
  20   /**
  21    * Widths of field in the full timestamp
  22    */
  23   enum {
  24     kBCWidth      = 12,
  25     kOrbitWidth   = 24,
  26     kPeriodWidth  = 28
  27   };
  28   /** Full time stamp */
  29   ULong64_t fFull;
  30   /** LDC creation time - standard time_t */
  31   UInt_t    fTime;
  32   /** Mask of detectors present in the event - FMD is 0x1000 */
  33   UInt_t    fDetectors;
  34   /** Type of event - 7 is physics */
  35   UShort_t  fType;
  36   /** GUID */
  37   UChar_t   fGUID[42];
  38   /**
  39    * Create a branch in a tree
  40    *
  41    * @param tree Tree to create the branch in
  42    */
  43   void CreateBranch(TTree* tree)
  44   {
  45     tree->Branch("event", &(this->fFull),
  46                  "full/l:time/i:detector:type/s:guid/C");
  47   }
  48   /**
  49    * Set the address of a branch for reading back objects from the tree
  50    *
  51    * @param tree Tree to read back from
  52    */
  53   void ReadBranch(TTree* tree)
  54   {
  55     tree->SetBranchAddress("event", &(this->fFull));
  56   }
  57   /**
  58    * Utility function to encode the full time stamp from components.
  59    *
  60    * @param period Period counter (overflow of orbit counter)
  61    * @param orbit  Orbit counter (period of 88us)
  62    * @param bc     Bunch crossing number (period of 25ns)
  63    *
  64    * @return Encoded full time stamp
  65    */
  66   static ULong64_t EncodeFull(ULong64_t period,
  67                               ULong64_t orbit,
  68                               ULong64_t bc)
  69   {
  70     const ULong64_t bcMask     = (1 << kBCWidth) - 1;
  71     const ULong64_t orbitMask  = (1 << kOrbitWidth) - 1;
  72     const ULong64_t periodMask = (1 << kPeriodWidth) - 1;
  73     ULong64_t ret = ((bc     & bcMask)                  |
  74                      (orbit  & orbitMask)  <<  kBCWidth |
  75                      (period & periodMask) << (kBCWidth+kOrbitWidth));
  76     return ret;
  77   }
  78   /**
  79    * Fill information from ESD into this data structure.
  80    *
  81    * @param esd  Event
  82    * @param dets List of active detectors in this event.
  83    * @param guid Current file GUID
  84    */
  85   void Fill(AliESDEvent* esd, UInt_t dets, const TString& guid);
  86 };
  87
  88 #ifndef NO_TASK
  89 # ifndef __CINT__
  90 #  include <AliESDEvent.h>
  91 # endif
  92 inline void EventTimeData::Fill(AliESDEvent* esd, UInt_t dets,
  93                                 const TString& guid)
  94 {
  95   ULong64_t period = esd->GetPeriodNumber();
  96   ULong64_t orbit  = esd->GetOrbitNumber();
  97   ULong64_t bc     = esd->GetBunchCrossNumber();
  98   fType            = esd->GetEventType();
  99   fTime            = esd->GetTimeStamp();//LDC time
 100   fDetectors       = dets; // esd->GetDAQDetectorPattern();
 101   fFull            = EncodeFull(period, orbit, bc);
 102   Int_t i = 0;
 103   for (i = 0; i < guid.Length() && i < 42; i++) fGUID[i] = guid[i];
 104   for (; i < 42; i++) fGUID[i] = '\0';
 105   fGUID[41] = '\0';
 106 }
 107 # else
 108 inline void EventTimeData::Fill(AliESDEvent*, UInt_t, const TString&)
 109 {
 110   Warning("Fill", "Calling empty method - shouldn't happen");
 111 }
 112 #endif
 113
 114 #ifndef NO_TASK
 115 # ifndef __CINT__
 116 #  include <AliAnalysisManager.h>
 117 #  include <AliVEventHandler.h>
 118 #  include <AliESDEvent.h>
 119 #  include <TTree.h>
 120 #  include <TH2.h>
 121 #  include <TList.h>
 122 #  include <AliCDBManager.h>
 123 #  include <AliCDBEntry.h>
 124 #  include <AliTriggerConfiguration.h>
 125 #  include <AliTriggerClass.h>
 126 #  include <AliTriggerCluster.h>
 127 #  include <AliDAQ.h>
 128 #  include <TObjArray.h>
 129 #  include <TDirectory.h>
 130 #  include <TUrl.h>
 131 # else
 132 class AliAnalysisManager;
 133 class TTree;
 134 class AliTimeStamp;
 135 class TList;
 136 class TH2;
 137 # endif
 138 # include <AliAnalysisTaskSE.h>
 139 # include <vector>
 140
 141 /**
 142  * A task to record the unique timestamp of each event.
 143  *
 144  * @par Input: ESD
 145  * @par Output:  A tree with a single branch
 146  * @ingroup pwglf_forward_eventtime
 147  */
 148 class EventTimeTask : public AliAnalysisTaskSE
 149 {
 150 public:
 151   enum {
 152     kListSlot = 1,
 153     kTreeSlot = 2
 154   };
 155   /**
 156    * Default CTOR - for I/O only.
 157    */
 158   EventTimeTask()
 159     : AliAnalysisTaskSE(),
 160       fTree(0),
 161       fHistograms(0),
 162       fDetVsType(0)
 163   {
 164   }
 165   /**
 166    * User constructor
 167    *
 168    * @param name Name of task
 169    */
 170   EventTimeTask(const char* name)
 171     : AliAnalysisTaskSE(name),
 172       fTree(0),
 173       fHistograms(0),
 174       fDetVsType(0),
 175       fGUID("")
 176   {
 177     DefineOutput(kListSlot, TList::Class());
 178     DefineOutput(kTreeSlot, TTree::Class());
 179     fBranchNames = "ESD:AliESDRun.,AliESDHeader.";
 180   }
 181   /**
 182    * Create user output objects.
 183    *
 184    * Called on each slave at start of job
 185    */
 186   void UserCreateOutputObjects()
 187   {
 188     Printf("Creating tree and histogram");
 189     fGUID = "";
 190     fHistograms = new TList();
 191     fHistograms->SetOwner();
 192     fHistograms->SetName("L");
 193
 194     fDetVsType = new TH2D("detVsType", "Detector vs type",
 195                      16, -.5, 15.5, 31, -.5, 30.5);
 196     fDetVsType->SetXTitle("Type");
 197     fDetVsType->SetYTitle("Detector");
 198     fDetVsType->SetDirectory(0);
 199     fHistograms->Add(fDetVsType);
 200     Printf("Histogram (%d,%f,%f)x(%d,%f,%f)",
 201            fDetVsType->GetXaxis()->GetNbins(),
 202            fDetVsType->GetXaxis()->GetXmin(),
 203            fDetVsType->GetXaxis()->GetXmax(),
 204            fDetVsType->GetYaxis()->GetNbins(),
 205            fDetVsType->GetYaxis()->GetXmin(),
 206            fDetVsType->GetYaxis()->GetXmax());
 207
 208     // TDirectory* savdir = gDirectory;
 209     // Printf("Opening file at slot %d", kTreeSlot);
 210     // OpenFile(kTreeSlot);
 211     Printf("Make tree and disassociate from file");
 212     fTree = new TTree("T", "T");
 213     fTree->SetDirectory(0);
 214     Printf("Create branch");
 215     fData.CreateBranch(fTree);
 216     // savdir->cd();
 217
 218
 219     PostData(kListSlot, fHistograms);
 220     PostData(kTreeSlot, fTree);
 221   }
 222   /**
 223    * Analyse a single event
 224    *
 225    */
 226   void UserExec(Option_t*)
 227   {
 228     static Bool_t first = true;
 229     LoadBranches();
 230
 231     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 232     if (!esd) return;
 233     if (!esd->GetHeader()) return;
 234
 235     if (first) {
 236       LoadTriggerConfig(esd->GetRunNumber());
 237       first = false;
 238     }
 239     ULong64_t mask = esd->GetTriggerMask();
 240     UInt_t    dets = 0;
 241     for (UShort_t i = 0; i < fDets.size(); i++) {
 242       if ((1 << i) & mask) dets |= fDets[i];
 243     }
 244     // Printf("Event mask 0x%016llx -> 0x%08x", mask, dets);
 245     fData.Fill(esd, dets, fGUID);
 246     fTree->Fill();
 247
 248     UInt_t type      = esd->GetEventType();
 249     UInt_t detectors = dets;
 250     for (UInt_t i = 0; i < 31; i++) {
 251       if ((1 << i) & detectors) fDetVsType->Fill(type, i);
 252     }
 253
 254     PostData(kListSlot, fHistograms);
 255     PostData(kTreeSlot, fTree);
 256   }
 257   Bool_t UserNotify()
 258   {
 259     fGUID = "";
 260     AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
 261     AliVEventHandler*   inp = mgr->GetInputEventHandler();
 262     if (!inp) {
 263       Warning("UserNotify", "No input handler");
 264       return true;
 265     }
 266     TTree* tree = inp->GetTree();
 267     if (!tree) {
 268       Warning("UserNotify", "No input tree");
 269       return true;
 270     }
 271     TFile* file = tree->GetCurrentFile();
 272     if (!file) {
 273       Warning("UserNotify", "No current file for tree");
 274       return true;
 275     }
 276     const TUrl* url = file->GetEndpointUrl();
 277     if (!url) {
 278       Warning("UserNotify", "No end point for file");
 279       return false;
 280     }
 281     fGUID = gSystem->BaseName(url->GetFile());
 282     Printf("Got GUID=%s from %s", fGUID.Data(), url->GetUrl());
 283
 284     return true;
 285   }
 286   void LoadTriggerConfig(Int_t runNo)
 287   {
 288     Printf("Loading trigger configuration for run %d",  runNo);
 289     // --- Connect to CDB --------------------------------------------
 290     AliCDBManager* cdb = AliCDBManager::Instance();
 291     cdb->SetDefaultStorageFromRun(runNo);
 292     cdb->SetRun(runNo);
 293
 294     // --- Get entry -------------------------------------------------
 295     AliCDBEntry* entry = cdb->Get("GRP/CTP/Config");
 296     if (!entry || !entry->GetObject()) {
 297       Warning("LoadTriggerConfig", "Couldn't get trigger configuration");
 298       return;
 299     }
 300     AliTriggerConfiguration* config =
 301       static_cast<AliTriggerConfiguration*>(entry->GetObject());
 302
 303     // --- Get the classes, and resize cache -------------------------
 304     const TObjArray& clss = config->GetClasses();
 305     fDets.resize(clss.GetEntries());
 306
 307     // --- Loop over configurations ----------------------------------
 308     TIter next(&clss);
 309     AliTriggerClass* cls = 0;
 310     while ((cls = static_cast<AliTriggerClass*>(next()))) {
 311       Int_t              mask    = cls->GetMask();
 312       AliTriggerCluster* cluster = cls->GetCluster();
 313       if (!cluster) {
 314         Warning("LoadTriggerConfig",
 315                 "Failed to get trigger cluster for %s", cls->GetName());
 316         continue;
 317       }
 318       TString names = cluster->GetDetectorsInCluster();
 319       if (names.IsNull()) {
 320         Warning("LoadTriggerConfig", "No detectors for cluster %s",
 321                 cls->GetName());
 322         continue;
 323       }
 324       UInt_t dets = AliDAQ::DetectorPattern(names);
 325       UInt_t id   = UInt_t(TMath::Log2(mask));
 326       fDets[id]   = dets;
 327       Printf("Trigger mask 0x%08x (%3d): 0x%08x (%s)",
 328              mask, id, dets, names.Data());
 329     }
 330     // cdb->Destroy();
 331   }
 332   /**
 333    * Register with manager and connect output containers
 334    *
 335    * @param mgr Manager
 336    */
 337   void Connect(AliAnalysisManager* mgr)
 338   {
 339     mgr->AddTask(this);
 340     mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
 341     AliAnalysisDataContainer* contTree =
 342       mgr->CreateContainer("T", TTree::Class(),
 343                            AliAnalysisManager::kOutputContainer,
 344                            "time.root");
 345     AliAnalysisDataContainer* contList =
 346       mgr->CreateContainer("L", TList::Class(),
 347                            AliAnalysisManager::kOutputContainer,
 348                            "hist.root");
 349     mgr->ConnectOutput(this, kTreeSlot, contTree);
 350     mgr->ConnectOutput(this, kListSlot, contList);
 351   }
 352   /**
 353    * Create an instance of this task, and register it and it's
 354    * outputs.
 355    */
 356   static void Create()
 357   {
 358     EventTimeTask* task = new EventTimeTask("eventTime");
 359     task->Connect(AliAnalysisManager::GetAnalysisManager());
 360   }
 361   TTree*              fTree;       // Our tree
 362   EventTimeData       fData;       // Our data
 363   TList*              fHistograms; // List
 364   TH2D*               fDetVsType;  // Histogram
 365   std::vector<UInt_t> fDets;       // Per-trigger-bit detector mask
 366   TString             fGUID;
 367
 368   ClassDef(EventTimeTask,3);
 369 };
 370 #endif // NO_TASK
 371
 372 #ifndef NO_MAP
 373 #include <utility>
 374 #include <map>
 375
 376 typedef std::pair<ULong64_t,ULong64_t> EventTimeMapPair;
 377
 378 /**
 379  * A map of event time-stamp to distance to previous event
 380  *
 381  * @ingroup pwglf_forward_eventtime
 382  */
 383 struct EventTimeMap : public TObject
 384 {
 385   /** Map type */
 386   typedef std::map<ULong64_t,ULong64_t> Map;
 387   /** Key type */
 388   typedef Map::key_type Key;
 389   /** Mapped value type */
 390   typedef Map::mapped_type Value;
 391   /** Element type */
 392   typedef Map::value_type Pair;
 393   /** Iterator type */
 394   typedef Map::iterator   Iterator;
 395   /** Constant iterator type */
 396   typedef Map::const_iterator ConstIterator;
 397   /**
 398    * Constructor
 399    */
 400   EventTimeMap() : TObject(), fMap() {}
 401   /**
 402    * Destructor
 403    */
 404   virtual ~EventTimeMap() {}
 405   /**
 406    * Get name of this object - always the same
 407    *
 408    * @return The string "eventTimeMap"
 409    */
 410   const char* GetName() const { return "eventTimeMap"; }
 411   /**
 412    * Element access.  If the key @a k doesn't already exist, it is
 413    * created
 414    *
 415    * @param k Key
 416    *
 417    * @return A pair of key and value
 418    */
 419   Value& operator[](const Key& k)
 420   {
 421     return fMap[k];
 422   }
 423   /**
 424    * Find the element whos key is @a k
 425    *
 426    * @param k Key to look for
 427    *
 428    * @return Iterator pointing to element, or end of container if not found
 429    */
 430   Iterator Find(const Key& k)
 431   {
 432     return fMap.find(k);
 433   }
 434   /**
 435    * Find the element whos key is @a k
 436    *
 437    * @param k Key to look for
 438    *
 439    * @return Iterator pointing to element, or end of container if not found
 440    */
 441   ConstIterator Find(const Key& k) const
 442   {
 443     return fMap.find(k);
 444   }
 445   /**
 446    * Get forward iterator pointing beginning of the container
 447    *
 448    * @return Iterator to start of container
 449    */
 450   Iterator Begin()
 451   {
 452     return fMap.begin();
 453   }
 454   /**
 455    * Get forward iterator pointing beginning of the container
 456    *
 457    * @return Iterator to start of container
 458    */
 459   ConstIterator Begin() const
 460   {
 461     return fMap.begin();
 462   }
 463   /**
 464    * Get forward iterator pointing just beyond the end of the container
 465    *
 466    * @return Iterator just beyond container
 467    */
 468   Iterator End()
 469   {
 470     return fMap.end();
 471   }
 472   /**
 473    * Get forward iterator pointing just beyond the end of the container
 474    *
 475    * @return Iterator just beyond container
 476    */
 477   ConstIterator End() const
 478   {
 479     return fMap.end();
 480   }
 481   enum {
 482     kInvalidTime = 0xFFFFFFFFFFFFFFFF
 483   };
 484   /**
 485    * Get the time difference to previous event from a event with a
 486    * given time stamp.
 487    *
 488    * @param timestamp Time stamp of the event
 489    *
 490    * @return time difference or kInvalidTime if not found
 491    */
 492   Value Get(const Key& timestamp) const
 493   {
 494     ConstIterator i = Find(timestamp);
 495     if (i == End()) return kInvalidTime;
 496     return i->second;
 497   }
 498   /**
 499    * Set the time to previous event for a given event time stamp
 500    *
 501    * @param timestamp Event time stamp
 502    * @param diff      Time to previous event
 503    */
 504   void Set(const Key& timestamp, const Value& diff)
 505   {
 506     this->operator[](timestamp) = diff;
 507   }
 508   ULong64_t Size() const { return fMap.size(); }
 509   /** Our embeded map */
 510   Map fMap;
 511   ClassDef(EventTimeMap,1)
 512 };
 513 #endif // NO_MAP
 514
 515 #ifndef NO_SORTER
 516 # ifndef __CINT__
 517 #  include <TFile.h>
 518 #  include <TArrayL64.h>
 519 #  include <TMath.h>
 520 #  include <TParameter.h>
 521 #  include <TCanvas.h>
 522 #  include <TH1.h>
 523 # else
 524 class TFile;
 525 class TTree;
 526 class TCanvas;
 527 # endif
 528
 529 /**
 530  * A class to sort the tree and generate our timestamp->dT map.
 531  *
 532  * @ingroup pwglf_forward_eventtime
 533  */
 534 struct EventTimeSorter
 535 {
 536   TTree*        fTree;
 537   EventTimeData fData;
 538   ULong64_t     fMin;
 539   ULong64_t     fMax;
 540
 541   /**
 542    * Constructor
 543    */
 544   EventTimeSorter() : fTree(0), fData(), fMin(0xFFFFFFFFFFFFFFFF), fMax(0) {}
 545   /**
 546    * Progress meter
 547    *
 548    * @param cur    Current entry
 549    * @param total  Total number of entries
 550    */
 551   void Progress(Long64_t cur, Long64_t total) const
 552   {
 553     static UShort_t old = 0;
 554     if (cur == 0) old = 0;
 555     UShort_t now = UShort_t(100 * (cur + 1 == total ? 1 :
 556                                    Double_t(cur)/total));
 557     if (now > old)
 558       std::cout << "\rLooping over " << total << " entries: ... "
 559                 << now << '%' << std::flush;
 560     if (now >= 100) std::cout << std::endl;
 561     old = now;
 562   }
 563   /**
 564    * Connect a tree
 565    *
 566    * @param inputName
 567    * @param treeName
 568    *
 569    * @return
 570    */
 571   Bool_t OpenInput(const char* inputName, const char* treeName)
 572   {
 573     CloseInput();
 574
 575     // --- Get input -------------------------------------------------
 576     TChain* chain = new TChain(treeName);
 577     chain->Add(inputName);
 578     if (chain->GetListOfFiles()->GetEntries() < 1) {
 579       Error("Run", "Failed to add \"%s\" to chain", inputName);
 580       return false;
 581     }
 582     fTree = chain;
 583
 584     // --- Set branch address ---------------------------------------
 585     fData.ReadBranch(fTree);
 586
 587     return true;
 588   }
 589   /**
 590    * Disconnect tree
 591    */
 592   void CloseInput()
 593   {
 594     if (!fTree) return;
 595     TFile* file = fTree->GetCurrentFile();
 596     if (file) file->Close();
 597     fTree = 0;
 598   }
 599   /**
 600    * Run the sorter
 601    *
 602    * @param inputName   Input file name
 603    * @param outputName  Output file name
 604    * @param treeName    Tree name
 605    *
 606    * @return true on success
 607    */
 608   Bool_t Run(const char* inputName, const char* outputName,
 609              const char* treeName="T")
 610   {
 611     if (!OpenInput(inputName, treeName)) return false;
 612
 613     // --- Loop over the data ----------------------------------------
 614     Long64_t   nEnt    = fTree->GetEntries();
 615     Long64_t   n       = 0;
 616     ULong64_t* values  = new ULong64_t[nEnt];
 617     UInt_t*    dets    = new UInt_t[nEnt];
 618     UShort_t*  types   = new UShort_t[nEnt];
 619     UInt_t*    times   = new UInt_t[nEnt];
 620     for (Long64_t iEnt = 0; iEnt < nEnt; iEnt++) {
 621       Progress(iEnt, nEnt);
 622
 623       fTree->GetEntry(iEnt);
 624
 625       if (!(fData.fDetectors & 0x1000))
 626         // No FMD read-out
 627         continue;
 628       if (fData.fType != 7) {
 629         // Ignore all but physics events
 630         Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType);
 631         // continue;
 632       }
 633
 634       // Store values
 635       Int_t j = n;
 636       values[j] = fData.fFull;
 637       dets[j]   = fData.fDetectors;
 638       types[j]  = fData.fType;
 639       times[j]  = fData.fTime;
 640       n++;
 641     }
 642
 643     // --- Now sort all values ---------------------------------------
 644     TArrayL64 index(n);
 645     TMath::Sort(n, values, index.fArray, false);
 646
 647     // Maps the unique time to the distance to previous event
 648     EventTimeMap  eventTimeMap;
 649     ULong64_t     key   = values[index[0]];
 650     ULong64_t     prev  = 0;
 651     ULong64_t     dt    = key-prev;
 652     ULong64_t     min   = 0xFFFFFFFFFFFFFFFF;
 653     ULong64_t     max   = 0;
 654     eventTimeMap[key]   = dt;
 655     for (Long64_t i = 1; i < n; i++) {
 656       Long64_t j        = index[i];
 657       Long64_t k        = index[i-1];
 658       key               = values[j];
 659       prev              = values[k];
 660       dt                = key - prev;
 661       if (dt <= 0) {
 662         Printf("0x%016llx==0x%016llx -> dt=%10llu [%10lld %10lld] "
 663                "(det: 0x%08x 0x%08x  type: 0x%04x 0x%04x  time: %08d %08d)",
 664                key, prev, dt, j, k, dets[j], dets[k],
 665                types[j], types[k], times[j], times[k]);
 666         // continue;
 667       }
 668       eventTimeMap[key] = dt;
 669       min               = TMath::Min(dt, min);
 670       max               = TMath::Max(dt, max);
 671     }
 672     std::cout << "Range is [" << min << "," << max << ']' << std::endl;
 673
 674     TFile* output = TFile::Open(outputName, "RECREATE");
 675     if (!output) {
 676       Error("Run", "Failed to create output file \"%s\"", outputName);
 677       return false;
 678     }
 679     fMin = min;
 680     fMax = max;
 681     eventTimeMap.Write();
 682     output->Write();
 683     output->Close();
 684
 685     delete [] values;
 686
 687     CloseInput();
 688
 689     return true;
 690   }
 691   Bool_t Test(const char* inputName, const char* outputName,
 692               const char* treeName="T")
 693   {
 694     if (!OpenInput(inputName, treeName)) return false;
 695
 696     // --- Get our map -----------------------------------------------
 697     TFile* output = TFile::Open(outputName, "UPDATE");
 698     if (!output) {
 699       Error("Test", "Failed to open \"%s\"", outputName);
 700       return false;
 701     }
 702
 703     EventTimeMap* eventTimeMap =
 704       static_cast<EventTimeMap*>(output->Get("eventTimeMap"));
 705     if (!eventTimeMap) {
 706       Error("Test", "Failed to get \"%s\" from \"%s\"",
 707             "eventTimeMap", outputName);
 708       return false;
 709     }
 710
 711     // --- Histograms --------------------------------------------------
 712     ULong64_t mmin = TMath::Min(25*fMin, 900000ULL);
 713     ULong64_t mmax = TMath::Min(25*fMax, 110000000ULL);
 714     TH1*      diff = new TH1D("diff", "Time-to-last-event (10#mus bins w/FMD)",
 715                               (mmax-mmin)/10000, mmin, mmax);
 716     diff->SetStats(0);
 717     diff->SetXTitle("#Deltat [ns]");
 718     diff->SetFillColor(kRed+1);
 719     diff->SetFillStyle(3001);
 720     diff->SetDirectory(0);
 721
 722
 723     TH1* ldiff = new TH1D("ldiff", "log(#Deltat) - Events w/FMD",
 724                           300, 0, 15);
 725     ldiff->SetStats(0);
 726     ldiff->SetXTitle("log_{10}(#Deltat) [ns]");
 727     ldiff->SetFillColor(kGreen+1);
 728     ldiff->SetFillStyle(3001);
 729     ldiff->SetDirectory(0);
 730
 731     // --- Loop over the data ----------------------------------------
 732     Long64_t   nEnt    = fTree->GetEntries();
 733     Long64_t   nZero   = 0;
 734     Long64_t   nMiss   = 0;
 735     Long64_t   nGood   = 0;
 736     Long64_t   nNoFMD  = 0;
 737     for (Long64_t iEnt = 0; iEnt < /*10*/ nEnt; iEnt++) {
 738       Progress(iEnt, nEnt);
 739       fTree->GetEntry(iEnt);
 740
 741       if (!(fData.fDetectors & 0x1000)) {
 742         // No FMD read-out
 743         nNoFMD++;
 744         continue;
 745       }
 746       if (fData.fType != 7) {
 747         // Ignore all but physics events
 748         Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType);
 749         // continue;
 750       }
 751
 752       // Look-up the timestamp
 753       ULong64_t value = fData.fFull;
 754       ULong64_t dT    = eventTimeMap->Get(value);
 755       if (dT == EventTimeMap::kInvalidTime) {
 756         Warning("Test", "Value %llu not found", value);
 757         ldiff->Fill(1);
 758         nMiss++;
 759         continue;
 760       }
 761       if (dT <= 0) {
 762 #if 0
 763         Warning("Test", "Impossible dt=%llu found for %llu (0x%0x %2d)",
 764                 dT, value, fData.fDetectors, fData.fType);
 765 #endif
 766         ldiff->Fill(1);
 767         nZero++;
 768         continue;
 769       }
 770       nGood++;
 771       Double_t  dt    = 25.*dT;
 772       diff->Fill(dt);
 773       Double_t  logDt = TMath::Log10(dt);
 774       ldiff->Fill(logDt);
 775     }
 776     CloseInput();
 777     Printf("missing: %llu, bad: %llu, good: %llu, no FMD: %llu, all: %llu",
 778            nMiss,nZero,nGood,nNoFMD,nEnt);
 779
 780     // --- Draw results ----------------------------------------------
 781     TCanvas* c = new TCanvas("c", "C");
 782     c->SetTopMargin(0.01);
 783     c->SetRightMargin(0.01);
 784     c->Divide(2,1); // ,0,0);
 785     TVirtualPad* p = c->cd(2);
 786     // p->SetRightMargin(0.10);
 787     p->SetLogy();
 788     ldiff->DrawCopy();
 789
 790     p = c->cd(1);
 791     p->SetLogy();
 792     p->SetLogx();
 793     diff->DrawCopy();
 794
 795     c->Print("dt.png");
 796     c->Print("dt.C");
 797
 798     // --- Write histogram to file -----------------------------------
 799     output->cd();
 800     diff->Write();
 801     output->Write();
 802     output->Close();
 803
 804     return true;
 805   }
 806 };
 807 #endif // NO_SORTER
 808 #ifdef __MAKECINT__
 809 # ifndef NO_MAP
 810 #  pragma link C++ class std::pair<ULong64_t,ULong64_t>+;
 811 #  pragma link C++ class std::map<ULong64_t,ULong64_t>+;
 812 # endif
 813 # ifndef NO_TASK
 814 #  pragma link C++ class std::vector<UInt_t>+;
 815 # endif
 816 #endif
 817
 818 // EOF
 819
 820