PWGLF/FORWARD/analysis2/trains/EventTimeTask.C

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