FMD/AliFMDInput.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16 /** @file    AliFMDInput.cxx
  17     @author  Christian Holm Christensen <cholm@nbi.dk>
  18     @date    Mon Mar 27 12:42:40 2006
  19     @brief   FMD utility classes for reading FMD data
  20 */
  21 //___________________________________________________________________
  22 //
  23 // The classes defined here, are utility classes for reading in data
  24 // for the FMD.  They are  put in a seperate library to not polute the
  25 // normal libraries.  The classes are intended to be used as base
  26 // classes for customized class that do some sort of analysis on the
  27 // various types of data produced by the FMD.
  28 //
  29 // Latest changes by Christian Holm Christensen
  30 //
  31 #include "AliFMDInput.h"        // ALIFMDHIT_H
  32 #include "AliFMDDebug.h"        // ALIFMDDEBUG_H ALILOG_H
  33 #include "AliLoader.h"          // ALILOADER_H
  34 #include "AliRunLoader.h"       // ALIRUNLOADER_H
  35 #include "AliRun.h"             // ALIRUN_H
  36 #include "AliStack.h"           // ALISTACK_H
  37 #include "AliRawReaderFile.h"   // ALIRAWREADERFILE_H
  38 #include "AliRawReaderRoot.h"   // ALIRAWREADERROOT_H
  39 #include "AliRawReaderDate.h"   // ALIRAWREADERDATE_H
  40 #include "AliRawEventHeaderBase.h"
  41 #include "AliFMD.h"             // ALIFMD_H
  42 #include "AliFMDHit.h"          // ALIFMDHIT_H
  43 #include "AliFMDDigit.h"        // ALIFMDDigit_H
  44 #include "AliFMDSDigit.h"       // ALIFMDDigit_H
  45 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
  46 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
  47 #include "AliFMDGeometry.h"
  48 #include <AliESD.h>
  49 #include <AliESDFMD.h>
  50 #include <AliESDEvent.h>
  51 #include <AliCDBManager.h>
  52 #include <AliCDBEntry.h>
  53 #include <AliAlignObjParams.h>
  54 #include <AliTrackReference.h>
  55 #include <TTree.h>              // ROOT_TTree
  56 #include <TChain.h>             // ROOT_TChain
  57 #include <TParticle.h>          // ROOT_TParticle
  58 #include <TString.h>            // ROOT_TString
  59 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
  60 #include <TMath.h>              // ROOT_TMath
  61 #include <TGeoManager.h>        // ROOT_TGeoManager
  62 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
  63 #include <Riostream.h>          // ROOT_Riostream
  64 #include <TFile.h>              // ROOT_TFile
  65 #include <TStreamerInfo.h>
  66 #include <TArrayF.h>
  67
  68 //____________________________________________________________________
  69 ClassImp(AliFMDInput)
  70 #if 0
  71   ; // This is here to keep Emacs for indenting the next line
  72 #endif
  73
  74 //____________________________________________________________________
  75 const AliFMDInput::ETrees AliFMDInput::fgkAllLoads[] = { kHits,
  76                                                          kKinematics,
  77                                                          kDigits,
  78                                                          kSDigits,
  79                                                          kHeader,
  80                                                          kRecPoints,
  81                                                          kESD,
  82                                                          kRaw,
  83                                                          kGeometry,
  84                                                          kTrackRefs,
  85                                                          kRawCalib,
  86                                                          kUser };
  87
  88 //____________________________________________________________________
  89 AliFMDInput::AliFMDInput()
  90   : TNamed("AliFMDInput", "Input handler for various FMD data"),
  91     fGAliceFile(""),
  92     fLoader(0),
  93     fRun(0),
  94     fStack(0),
  95     fFMDLoader(0),
  96     fReader(0),
  97     fFMDReader(0),
  98     fFMD(0),
  99     fESD(0),
 100     fESDEvent(0),
 101     fTreeE(0),
 102     fTreeH(0),
 103     fTreeTR(0),
 104     fTreeD(0),
 105     fTreeS(0),
 106     fTreeR(0),
 107     fTreeA(0),
 108     fChainE(0),
 109     fArrayE(0),
 110     fArrayH(0),
 111     fArrayTR(0),
 112     fArrayD(0),
 113     fArrayS(0),
 114     fArrayR(0),
 115     fArrayA(0),
 116     fHeader(0),
 117     fGeoManager(0),
 118     fTreeMask(0),
 119     fRawFile(""),
 120     fInputDir("."),
 121     fIsInit(kFALSE),
 122     fEventCount(0),
 123     fNEvents(-1)
 124 {
 125
 126   // Constructor of an FMD input object.  Specify what data to read in
 127   // using the AddLoad member function.   Sub-classes should at a
 128   // minimum overload the member function Event.   A full job can be
 129   // executed using the member function Run.
 130 }
 131
 132
 133
 134 //____________________________________________________________________
 135 AliFMDInput::AliFMDInput(const char* gAliceFile)
 136   : TNamed("AliFMDInput", "Input handler for various FMD data"),
 137     fGAliceFile(gAliceFile),
 138     fLoader(0),
 139     fRun(0),
 140     fStack(0),
 141     fFMDLoader(0),
 142     fReader(0),
 143     fFMDReader(0),
 144     fFMD(0),
 145     fESD(0),
 146     fESDEvent(0),
 147     fTreeE(0),
 148     fTreeH(0),
 149     fTreeTR(0),
 150     fTreeD(0),
 151     fTreeS(0),
 152     fTreeR(0),
 153     fTreeA(0),
 154     fChainE(0),
 155     fArrayE(0),
 156     fArrayH(0),
 157     fArrayTR(0),
 158     fArrayD(0),
 159     fArrayS(0),
 160     fArrayR(0),
 161     fArrayA(0),
 162     fHeader(0),
 163     fGeoManager(0),
 164     fTreeMask(0),
 165     fRawFile(""),
 166     fInputDir("."),
 167     fIsInit(kFALSE),
 168     fEventCount(0),
 169     fNEvents(-1)
 170 {
 171
 172   // Constructor of an FMD input object.  Specify what data to read in
 173   // using the AddLoad member function.   Sub-classes should at a
 174   // minimum overload the member function Event.   A full job can be
 175   // executed using the member function Run.
 176 }
 177
 178 //____________________________________________________________________
 179 void
 180 AliFMDInput::SetLoads(UInt_t mask)
 181 {
 182   for (UInt_t i = 0; i < sizeof(mask); i++) {
 183     if (!(mask & (1 << i))) continue;
 184     const ETrees *ptype = fgkAllLoads;
 185     do {
 186       ETrees type = *ptype;
 187       if (i != UInt_t(type)) continue;
 188       AddLoad(type);
 189       break;
 190     } while (*ptype++ != kUser);
 191   }
 192 }
 193
 194 //____________________________________________________________________
 195 void
 196 AliFMDInput::SetLoads(const char* what)
 197 {
 198   TString    l(what);
 199   TObjArray* ll = l.Tokenize(", ");
 200   TIter      next(ll);
 201   TObject*   os = 0;
 202   while ((os = next())) {
 203     ETrees type = ParseLoad(os->GetName());
 204     AddLoad(type);
 205   }
 206 }
 207
 208
 209 //____________________________________________________________________
 210 AliFMDInput::ETrees
 211 AliFMDInput::ParseLoad(const char* what)
 212 {
 213   TString opt(what);
 214   opt.ToLower();
 215   const ETrees* ptype = fgkAllLoads;
 216   do {
 217     ETrees  type = *ptype;
 218     if (opt.Contains(TreeName(type,true), TString::kIgnoreCase))
 219       return type;
 220   } while (*ptype++ != kUser);
 221   return kUser;
 222 }
 223 //____________________________________________________________________
 224 const char*
 225 AliFMDInput::LoadedString(Bool_t dataOnly) const
 226 {
 227   static TString ret;
 228   if (!ret.IsNull()) return ret.Data();
 229
 230   const ETrees* ptype = fgkAllLoads;
 231   do {
 232     ETrees type = *ptype;
 233     if (dataOnly &&
 234         (type == kKinematics ||
 235          type == kHeader     ||
 236          type == kGeometry   ||
 237          type == kTrackRefs)) continue;
 238     if (!IsLoaded(*ptype)) continue;
 239
 240     if (!ret.IsNull()) ret.Append(",");
 241     ret.Append(TreeName(type));
 242   } while (*ptype++ != kUser);
 243   return ret.Data();
 244 }
 245
 246 //____________________________________________________________________
 247 const char*
 248 AliFMDInput::TreeName(ETrees tree, Bool_t shortest)
 249 {
 250   if (shortest) {
 251     switch (tree) {
 252     case kHits:         return "hit";
 253     case kKinematics:   return "kin";
 254     case kDigits:       return "dig";
 255     case kSDigits:      return "sdig";
 256     case kHeader:       return "hea";
 257     case kRecPoints:    return "recp";
 258     case kESD:          return "esd";
 259     case kRaw:          return "raw";
 260     case kGeometry:     return "geo";
 261     case kTrackRefs:    return "trackr";
 262     case kRawCalib:     return "rawc";
 263     case kUser:         return "user";
 264     }
 265     return 0;
 266   }
 267   switch (tree) {
 268   case kHits:           return "Hits";
 269   case kKinematics:     return "Kinematics";
 270   case kDigits:         return "Digits";
 271   case kSDigits:        return "SDigits";
 272   case kHeader:         return "Header";
 273   case kRecPoints:      return "RecPoints";
 274   case kESD:            return "ESD";
 275   case kRaw:            return "Raw";
 276   case kGeometry:       return "Geometry";
 277   case kTrackRefs:      return "TrackRefs";
 278   case kRawCalib:       return "RawCalib";
 279   case kUser:           return "User";
 280   }
 281   return 0;
 282 }
 283
 284 //____________________________________________________________________
 285 Int_t
 286 AliFMDInput::NEvents() const
 287 {
 288   // Get number of events
 289   if (IsLoaded(kRaw) ||
 290       IsLoaded(kRawCalib)) return fReader->GetNumberOfEvents();
 291   if (fChainE) return fChainE->GetEntriesFast();
 292   if (fTreeE) return fTreeE->GetEntries();
 293   return -1;
 294 }
 295
 296 //____________________________________________________________________
 297 Bool_t
 298 AliFMDInput::Init()
 299 {
 300   // Initialize the object.  Get the needed loaders, and such.
 301
 302   // Check if we have been initialized
 303   if (fIsInit) {
 304     AliWarning("Already initialized");
 305     return fIsInit;
 306   }
 307   TString what;
 308   const ETrees* ptype = fgkAllLoads;
 309   do {
 310     ETrees type = *ptype;
 311     what.Append(Form("\n\t%-20s: %s", TreeName(type),
 312                      IsLoaded(type) ? "yes" : "no"));
 313   } while (*ptype++ != kUser);
 314
 315   Info("Init","Initialising w/mask 0x%04x%s", fTreeMask, what.Data());
 316   // Get the run
 317   if (IsLoaded(kDigits)     ||
 318       IsLoaded(kSDigits)    ||
 319       IsLoaded(kKinematics) ||
 320       IsLoaded(kTrackRefs)  ||
 321       IsLoaded(kHeader)) {
 322     if (!gSystem->FindFile(".:/", fGAliceFile)) {
 323       AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
 324     }
 325     else {
 326       fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
 327       if (!fLoader) {
 328         AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
 329         return kFALSE;
 330       }
 331       AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
 332
 333       if  (fLoader->LoadgAlice()) return kFALSE;
 334
 335       fRun = fLoader->GetAliRun();
 336
 337       // Get the FMD
 338       fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
 339       if (!fFMD) {
 340         AliError("Failed to get detector FMD from loader");
 341         return kFALSE;
 342       }
 343
 344       // Get the FMD loader
 345       fFMDLoader = fLoader->GetLoader("FMDLoader");
 346       if (!fFMDLoader) {
 347         AliError("Failed to get detector FMD loader from loader");
 348         return kFALSE;
 349       }
 350       if (fLoader->LoadHeader()) {
 351         AliError("Failed to get event header information from loader");
 352         return kFALSE;
 353       }
 354       fTreeE = fLoader->TreeE();
 355     }
 356   }
 357
 358   // Optionally, get the ESD files
 359   if (IsLoaded(kESD)) {
 360     fChainE = MakeChain("ESD", fInputDir, true);
 361     fESDEvent = new AliESDEvent();
 362     fESDEvent->ReadFromTree(fChainE);
 363     //    fChainE->SetBranchAddress("ESD", &fMainESD);
 364
 365   }
 366
 367   if (IsLoaded(kRaw) ||
 368       IsLoaded(kRawCalib)) {
 369     AliInfo("Getting FMD raw data digits");
 370     fArrayA = new TClonesArray("AliFMDDigit");
 371 #if 0
 372     if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
 373       fReader = new AliRawReaderRoot(fRawFile.Data());
 374     else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
 375       fReader = new AliRawReaderDate(fRawFile.Data());
 376     else
 377       fReader = new AliRawReaderFile(-1);
 378 #else
 379     if(!fRawFile.IsNull())
 380       fReader = AliRawReader::Create(fRawFile.Data());
 381     else
 382       fReader = new AliRawReaderFile(-1);
 383 #endif
 384     fFMDReader = new AliFMDRawReader(fReader, 0);
 385   }
 386
 387   // Optionally, get the geometry
 388   if (IsLoaded(kGeometry)) {
 389     TString fname;
 390     if (fRun) {
 391       fname = gSystem->DirName(fGAliceFile);
 392       fname.Append("/geometry.root");
 393     }
 394     if (!gSystem->AccessPathName(fname.Data()))
 395       fname = "";
 396     AliCDBManager* cdb   = AliCDBManager::Instance();
 397     if (!cdb->IsDefaultStorageSet()) {
 398       cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 399       cdb->SetRun(0);
 400     }
 401
 402     AliGeomManager::LoadGeometry(fname.IsNull() ? 0 : fname.Data());
 403
 404     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
 405     if (align) {
 406       AliInfo("Got alignment data from CDB");
 407       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
 408       if (!array) {
 409         AliWarning("Invalid align data from CDB");
 410       }
 411       else {
 412         Int_t nAlign = array->GetEntries();
 413         for (Int_t i = 0; i < nAlign; i++) {
 414           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
 415           if (!a->ApplyToGeometry()) {
 416             AliWarning(Form("Failed to apply alignment to %s",
 417                             a->GetSymName()));
 418           }
 419         }
 420       }
 421     }
 422     AliFMDGeometry* geom = AliFMDGeometry::Instance();
 423     geom->Init();
 424     geom->InitTransformations();
 425   }
 426
 427   fEventCount = 0;
 428   fIsInit = kTRUE;
 429   return fIsInit;
 430 }
 431
 432 //____________________________________________________________________
 433 Bool_t
 434 AliFMDInput::Begin(Int_t event)
 435 {
 436   // Called at the begining of each event.  Per default, it gets the
 437   // data trees and gets pointers to the output arrays.   Users can
 438   // overload this, but should call this member function in the
 439   // overloaded member function of the derived class.
 440
 441   // Check if we have been initialized
 442   if (!fIsInit) {
 443     AliError("Not initialized");
 444     return fIsInit;
 445   }
 446
 447   // Get the event
 448   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
 449
 450   // Possibly load global kinematics information
 451   if (IsLoaded(kKinematics)) {
 452     // AliInfo("Getting kinematics");
 453     if (fLoader->LoadKinematics("READ")) return kFALSE;
 454     fStack = fLoader->Stack();
 455   }
 456
 457   // Possibly load FMD Hit information
 458   if (IsLoaded(kHits)) {
 459     // AliInfo("Getting FMD hits");
 460     if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
 461     fTreeH = fFMDLoader->TreeH();
 462     if (!fArrayH) fArrayH = fFMD->Hits();
 463   }
 464
 465   // Possibly load FMD TrackReference information
 466   if (IsLoaded(kTrackRefs)) {
 467     // AliInfo("Getting FMD hits");
 468     if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
 469     fTreeTR = fLoader->TreeTR();
 470     if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
 471     fTreeTR->SetBranchAddress("TrackReferences",  &fArrayTR);
 472   }
 473
 474   // Possibly load heaedr information
 475   if (IsLoaded(kHeader)) {
 476     // AliInfo("Getting FMD hits");
 477     if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
 478     fHeader = fLoader->GetHeader();
 479   }
 480
 481   // Possibly load FMD Digit information
 482   if (IsLoaded(kDigits)) {
 483     // AliInfo("Getting FMD digits");
 484     if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
 485     fTreeD = fFMDLoader->TreeD();
 486     if (fTreeD) {
 487       if (!fArrayD) fArrayD = fFMD->Digits();
 488     }
 489     else {
 490       fArrayD = 0;
 491       AliWarning(Form("Failed to load FMD Digits"));
 492     }
 493   }
 494
 495   // Possibly load FMD Sdigit information
 496   if (IsLoaded(kSDigits)) {
 497     // AliInfo("Getting FMD summable digits");
 498     if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) {
 499       AliWarning("Failed to load SDigits!");
 500       return kFALSE;
 501     }
 502     fTreeS = fFMDLoader->TreeS();
 503     if (!fArrayS) fArrayS = fFMD->SDigits();
 504   }
 505
 506   // Possibly load FMD RecPoints information
 507   if (IsLoaded(kRecPoints)) {
 508     // AliInfo("Getting FMD reconstructed points");
 509     if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
 510     fTreeR = fFMDLoader->TreeR();
 511     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
 512     fTreeR->SetBranchAddress("FMD",  &fArrayR);
 513   }
 514
 515   // Possibly load FMD ESD information
 516   if (IsLoaded(kESD)) {
 517     // AliInfo("Getting FMD event summary data");
 518     Int_t read = fChainE->GetEntry(event);
 519     if (read <= 0) return kFALSE;
 520     fESD = fESDEvent->GetFMDData();
 521     if (!fESD) return kFALSE;
 522   }
 523
 524   // Possibly load FMD Digit information
 525   if (IsLoaded(kRaw) || IsLoaded(kRawCalib)) {
 526     Bool_t mon = fRawFile.Contains("mem://");
 527     // AliInfo("Getting FMD raw data digits");
 528     if (mon) std::cout << "Waiting for event ..." << std::flush;
 529     do {
 530       if (!fReader->NextEvent()) {
 531         if (mon) {
 532           gSystem->Sleep(3);
 533           continue;
 534         }
 535         return kFALSE;
 536       }
 537       UInt_t eventType = fReader->GetType();
 538       if(eventType == AliRawEventHeaderBase::kPhysicsEvent ||
 539          eventType == AliRawEventHeaderBase::kCalibrationEvent)
 540         break;
 541     } while (true);
 542     if (mon) std::cout << "got it" << std::endl;
 543     // AliFMDRawReader r(fReader, 0);
 544     fArrayA->Clear();
 545     fFMDReader->ReadAdcs(fArrayA);
 546     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
 547   }
 548   fEventCount++;
 549   return kTRUE;
 550 }
 551
 552
 553 //____________________________________________________________________
 554 Bool_t
 555 AliFMDInput::Event()
 556 {
 557   // Process one event.  The default implementation one or more of
 558   //
 559   //   -  ProcessHits       if the hits are loaded.
 560   //   -  ProcessDigits     if the digits are loaded.
 561   //   -  ProcessSDigits    if the sumbable digits are loaded.
 562   //   -  ProcessRecPoints  if the reconstructed points are loaded.
 563   //   -  ProcessESD        if the event summary data is loaded
 564   //
 565   if (IsLoaded(kHits))     if (!ProcessHits())      return kFALSE;
 566   if (IsLoaded(kTrackRefs))if (!ProcessTrackRefs()) return kFALSE;
 567   if (IsLoaded(kKinematics) &&
 568       IsLoaded(kHits))     if (!ProcessTracks())    return kFALSE;
 569   if (IsLoaded(kKinematics))if (!ProcessStack())     return kFALSE;
 570   if (IsLoaded(kSDigits))  if (!ProcessSDigits())   return kFALSE;
 571   if (IsLoaded(kDigits))   if (!ProcessDigits())    return kFALSE;
 572   if (IsLoaded(kRaw))      if (!ProcessRawDigits()) return kFALSE;
 573   if (IsLoaded(kRawCalib)) if (!ProcessRawCalibDigits())return kFALSE;
 574   if (IsLoaded(kRecPoints))if (!ProcessRecPoints()) return kFALSE;
 575   if (IsLoaded(kESD))      if (!ProcessESDs())      return kFALSE;
 576   if (IsLoaded(kUser))     if (!ProcessUsers())     return kFALSE;
 577
 578   return kTRUE;
 579 }
 580
 581 //____________________________________________________________________
 582 Bool_t
 583 AliFMDInput::ProcessHits()
 584 {
 585   // Read the hit tree, and pass each hit to the member function
 586   // ProcessHit.
 587   if (!fTreeH) {
 588     AliError("No hit tree defined");
 589     return kFALSE;
 590   }
 591   if (!fArrayH) {
 592     AliError("No hit array defined");
 593     return kFALSE;
 594   }
 595
 596   Int_t nTracks = fTreeH->GetEntries();
 597   for (Int_t i = 0; i < nTracks; i++) {
 598     Int_t hitRead  = fTreeH->GetEntry(i);
 599     if (hitRead <= 0) continue;
 600
 601     Int_t nHit = fArrayH->GetEntries();
 602     if (nHit <= 0) continue;
 603
 604     for (Int_t j = 0; j < nHit; j++) {
 605       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 606       if (!hit) continue;
 607
 608       TParticle* track = 0;
 609       if (IsLoaded(kKinematics) && fStack) {
 610         Int_t trackno = hit->Track();
 611         track = fStack->Particle(trackno);
 612       }
 613       if (!ProcessHit(hit, track)) return kFALSE;
 614     }
 615   }
 616   return kTRUE;
 617 }
 618 //____________________________________________________________________
 619 Bool_t
 620 AliFMDInput::ProcessTrackRefs()
 621 {
 622   // Read the reconstrcted points tree, and pass each reconstruction
 623   // object (AliFMDRecPoint) to either ProcessRecPoint.
 624   if (!fTreeTR) {
 625     AliError("No track reference tree defined");
 626     return kFALSE;
 627   }
 628   if (!fArrayTR) {
 629     AliError("No track reference array defined");
 630     return kFALSE;
 631   }
 632
 633   Int_t nEv = fTreeTR->GetEntries();
 634   for (Int_t i = 0; i < nEv; i++) {
 635     Int_t trRead  = fTreeTR->GetEntry(i);
 636     if (trRead <= 0) continue;
 637     Int_t nTrackRefs = fArrayTR->GetEntries();
 638     for (Int_t j = 0; j < nTrackRefs; j++) {
 639       AliTrackReference* trackRef =
 640         static_cast<AliTrackReference*>(fArrayTR->At(j));
 641       if (!trackRef) continue;
 642       // if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
 643       TParticle* track = 0;
 644       if (IsLoaded(kKinematics) && fStack) {
 645         Int_t trackno = trackRef->GetTrack();
 646         track = fStack->Particle(trackno);
 647       }
 648       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
 649     }
 650   }
 651   return kTRUE;
 652 }
 653 //____________________________________________________________________
 654 Bool_t
 655 AliFMDInput::ProcessTracks()
 656 {
 657   // Read the hit tree, and pass each hit to the member function
 658   // ProcessTrack.
 659   if (!fStack) {
 660     AliError("No track tree defined");
 661     return kFALSE;
 662   }
 663   if (!fTreeH) {
 664     AliError("No hit tree defined");
 665     return kFALSE;
 666   }
 667   if (!fArrayH) {
 668     AliError("No hit array defined");
 669     return kFALSE;
 670   }
 671
 672   // Int_t nTracks = fStack->GetNtrack();
 673   Int_t nTracks = fTreeH->GetEntries();
 674   for (Int_t i = 0; i < nTracks; i++) {
 675     Int_t      trackno = nTracks - i - 1;
 676     TParticle* track   = fStack->Particle(trackno);
 677     if (!track) continue;
 678
 679     // Get the hits for this track.
 680     Int_t hitRead  = fTreeH->GetEntry(i);
 681     Int_t nHit     = fArrayH->GetEntries();
 682     if (nHit == 0 || hitRead <= 0) {
 683       // Let user code see the track, even if there's no hits.
 684       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
 685       continue;
 686     }
 687
 688     // Loop over the hits corresponding to this track.
 689     for (Int_t j = 0; j < nHit; j++) {
 690       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 691       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
 692     }
 693   }
 694   return kTRUE;
 695 }
 696 //____________________________________________________________________
 697 Bool_t
 698 AliFMDInput::ProcessStack()
 699 {
 700   // Read the hit tree, and pass each hit to the member function
 701   // ProcessTrack.
 702   if (!fStack) {
 703     AliError("No track tree defined");
 704     return kFALSE;
 705   }
 706   Int_t nTracks = fStack->GetNtrack();
 707   for (Int_t i = 0; i < nTracks; i++) {
 708     Int_t      trackno = nTracks - i - 1;
 709     TParticle* track   = fStack->Particle(trackno);
 710     if (!track) continue;
 711
 712     if (!ProcessParticle(trackno, track)) return kFALSE;
 713   }
 714   return kTRUE;
 715 }
 716 //____________________________________________________________________
 717 Bool_t
 718 AliFMDInput::ProcessDigits()
 719 {
 720   // Read the digit tree, and pass each digit to the member function
 721   // ProcessDigit.
 722   if (!fTreeD) {
 723     AliError("No digit tree defined");
 724     return kFALSE;
 725   }
 726   if (!fArrayD) {
 727     AliError("No digit array defined");
 728     return kFALSE;
 729   }
 730
 731   Int_t nEv = fTreeD->GetEntries();
 732   for (Int_t i = 0; i < nEv; i++) {
 733     Int_t digitRead  = fTreeD->GetEntry(i);
 734     if (digitRead <= 0) continue;
 735     Int_t nDigit = fArrayD->GetEntries();
 736     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
 737     if (nDigit <= 0) continue;
 738     for (Int_t j = 0; j < nDigit; j++) {
 739       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 740       if (!digit) continue;
 741       if (!ProcessDigit(digit)) return kFALSE;
 742     }
 743   }
 744   return kTRUE;
 745 }
 746
 747 //____________________________________________________________________
 748 Bool_t
 749 AliFMDInput::ProcessSDigits()
 750 {
 751   // Read the summable digit tree, and pass each sumable digit to the
 752   // member function ProcessSdigit.
 753   if (!fTreeS) {
 754     AliWarning("No sdigit tree defined");
 755     return kTRUE; // Empty SDigits is fine
 756   }
 757   if (!fArrayS) {
 758     AliWarning("No sdigit array defined");
 759     return kTRUE; // Empty SDigits is fine
 760   }
 761
 762   Int_t nEv = fTreeS->GetEntries();
 763   for (Int_t i = 0; i < nEv; i++) {
 764     Int_t sdigitRead  = fTreeS->GetEntry(i);
 765     if (sdigitRead <= 0) {
 766       AliInfo(Form("Read nothing from tree"));
 767       continue;
 768     }
 769     Int_t nSdigit = fArrayS->GetEntriesFast();
 770     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
 771     AliInfo(Form("Got %5d digits for this event", nSdigit));
 772     if (nSdigit <= 0) continue;
 773     for (Int_t j = 0; j < nSdigit; j++) {
 774       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 775       if (!sdigit) continue;
 776       if (!ProcessSDigit(sdigit)) return kFALSE;
 777     }
 778   }
 779   return kTRUE;
 780 }
 781
 782 //____________________________________________________________________
 783 Bool_t
 784 AliFMDInput::ProcessRawDigits()
 785 {
 786   // Read the digit tree, and pass each digit to the member function
 787   // ProcessDigit.
 788   if (!fArrayA) {
 789     AliError("No raw digit array defined");
 790     return kFALSE;
 791   }
 792
 793   Int_t nDigit = fArrayA->GetEntries();
 794   if (nDigit <= 0) return kTRUE;
 795   for (Int_t j = 0; j < nDigit; j++) {
 796     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 797     if (!digit) continue;
 798     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 799       digit->Print();
 800     if (!ProcessRawDigit(digit)) return kFALSE;
 801   }
 802   return kTRUE;
 803 }
 804
 805 //____________________________________________________________________
 806 Bool_t
 807 AliFMDInput::ProcessRawCalibDigits()
 808 {
 809   // Read the digit tree, and pass each digit to the member function
 810   // ProcessDigit.
 811   if (!fArrayA) {
 812     AliError("No raw digit array defined");
 813     return kFALSE;
 814   }
 815
 816   Int_t nDigit = fArrayA->GetEntries();
 817   if (nDigit <= 0) return kTRUE;
 818   for (Int_t j = 0; j < nDigit; j++) {
 819     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 820     if (!digit) continue;
 821     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 822       digit->Print();
 823     if (!ProcessRawCalibDigit(digit)) return kFALSE;
 824   }
 825   return kTRUE;
 826 }
 827
 828 //____________________________________________________________________
 829 Bool_t
 830 AliFMDInput::ProcessRecPoints()
 831 {
 832   // Read the reconstrcted points tree, and pass each reconstruction
 833   // object (AliFMDRecPoint) to either ProcessRecPoint.
 834   if (!fTreeR) {
 835     AliError("No recpoint tree defined");
 836     return kFALSE;
 837   }
 838   if (!fArrayR) {
 839     AliError("No recpoints array defined");
 840     return kFALSE;
 841   }
 842
 843   Int_t nEv = fTreeR->GetEntries();
 844   for (Int_t i = 0; i < nEv; i++) {
 845     Int_t recRead  = fTreeR->GetEntry(i);
 846     if (recRead <= 0) continue;
 847     Int_t nRecPoint = fArrayR->GetEntries();
 848     for (Int_t j = 0; j < nRecPoint; j++) {
 849       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 850       if (!recPoint) continue;
 851       if (!ProcessRecPoint(recPoint)) return kFALSE;
 852     }
 853   }
 854   return kTRUE;
 855 }
 856
 857 //____________________________________________________________________
 858 Bool_t
 859 AliFMDInput::ProcessESDs()
 860 {
 861   // Process event summary data
 862   if (!fESD) return kFALSE;
 863   for (UShort_t det = 1; det <= 3; det++) {
 864     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 865     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 866       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 867       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 868       for (UShort_t sec = 0; sec < nsec; sec++) {
 869         for (UShort_t str = 0; str < nstr; str++) {
 870           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 871           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 872           if (!fESD->IsAngleCorrected())
 873             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 874           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 875         }
 876       }
 877     }
 878   }
 879   return kTRUE;
 880 }
 881
 882 //____________________________________________________________________
 883 Bool_t
 884 AliFMDInput::ProcessUsers()
 885 {
 886   // Process event summary data
 887   for (UShort_t det = 1; det <= 3; det++) {
 888     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 889     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 890       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 891       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 892       for (UShort_t sec = 0; sec < nsec; sec++) {
 893         for (UShort_t str = 0; str < nstr; str++) {
 894           Float_t v  = GetSignal(det,*rng,sec,str);
 895           if (!ProcessUser(det, *rng, sec, str, v)) continue;
 896         }
 897       }
 898     }
 899   }
 900   return kTRUE;
 901 }
 902
 903 //____________________________________________________________________
 904 Bool_t
 905 AliFMDInput::End()
 906 {
 907   // Called at the end of each event.  Per default, it unloads the
 908   // data trees and resets the pointers to the output arrays.   Users
 909   // can overload this, but should call this member function in the
 910   // overloaded member function of the derived class.
 911
 912   // Check if we have been initialized
 913   if (!fIsInit) {
 914     AliError("Not initialized");
 915     return fIsInit;
 916   }
 917   // Possibly unload global kinematics information
 918   if (IsLoaded(kKinematics)) {
 919     fLoader->UnloadKinematics();
 920     // fTreeK = 0;
 921     fStack = 0;
 922   }
 923   // Possibly unload FMD Hit information
 924   if (IsLoaded(kHits)) {
 925     fFMDLoader->UnloadHits();
 926     fTreeH = 0;
 927   }
 928   // Possibly unload FMD Digit information
 929   if (IsLoaded(kDigits)) {
 930     fFMDLoader->UnloadDigits();
 931     fTreeD = 0;
 932   }
 933   // Possibly unload FMD Sdigit information
 934   if (IsLoaded(kSDigits)) {
 935     fFMDLoader->UnloadSDigits();
 936     fTreeS = 0;
 937   }
 938   // Possibly unload FMD RecPoints information
 939   if (IsLoaded(kRecPoints)) {
 940     fFMDLoader->UnloadRecPoints();
 941     fTreeR = 0;
 942   }
 943   // AliInfo("Now out event");
 944   return kTRUE;
 945 }
 946
 947 //____________________________________________________________________
 948 Bool_t
 949 AliFMDInput::Run(UInt_t maxEvents)
 950 {
 951   // Run over all events and files references in galice.root
 952
 953   Bool_t retval;
 954   if (!(retval = Init())) return retval;
 955
 956   fNEvents = NEvents();
 957   if (fNEvents < 0)       fNEvents = maxEvents;
 958   else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
 959
 960   Int_t event = 0;
 961   for (; fNEvents < 0 || event < fNEvents; event++) {
 962     printf("\rEvent %8d/%8d ...", event, fNEvents);
 963     if (!(retval = Begin(event))) break;
 964     if (!(retval = Event())) break;
 965     if (!(retval = End())) break;
 966   }
 967   printf("Looped over %8d events\n", event+1);
 968   if (!retval) return retval;
 969   retval = Finish();
 970   return retval;
 971 }
 972
 973 //__________________________________________________________________
 974 TArrayF
 975 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
 976 {
 977   // Service function to define a logarithmic axis.
 978   // Parameters:
 979   //   n    Number of bins
 980   //   min  Minimum of axis
 981   //   max  Maximum of axis
 982   TArrayF bins(n+1);
 983   bins[0]      = min;
 984   if (n <= 20) {
 985     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
 986     return bins;
 987   }
 988   Float_t dp   = n / TMath::Log10(max / min);
 989   Float_t pmin = TMath::Log10(min);
 990   for (Int_t i = 1; i < n+1; i++) {
 991     Float_t p = pmin + i / dp;
 992     bins[i]   = TMath::Power(10, p);
 993   }
 994   return bins;
 995 }
 996
 997 //____________________________________________________________________
 998 void
 999 AliFMDInput::ScanDirectory(TSystemDirectory* dir,
1000                            const TString& olddir,
1001                            TChain* chain,
1002                            const char* pattern, bool recursive)
1003 {
1004   // Get list of files, and go back to old working directory
1005   TString oldDir(gSystem->WorkingDirectory());
1006   TList* files = dir->GetListOfFiles();
1007   gSystem->ChangeDirectory(oldDir);
1008
1009   // Sort list of files and check if we should add it
1010   if (!files) return;
1011   files->Sort();
1012   TIter next(files);
1013   TSystemFile* file = 0;
1014   while ((file = static_cast<TSystemFile*>(next()))) {
1015     TString name(file->GetName());
1016
1017     // Ignore special links
1018     if (name == "." || name == "..") continue;
1019
1020     // Check if this is a directory
1021     if (file->IsDirectory()) {
1022       if (recursive)
1023         ScanDirectory(static_cast<TSystemDirectory*>(file),
1024                       olddir, chain,
1025                       pattern,recursive);
1026       continue;
1027     }
1028
1029     // If this is not a root file, ignore
1030     if (!name.EndsWith(".root")) continue;
1031
1032     // If this file does not contain the pattern, ignore
1033     if (!name.Contains(pattern)) continue;
1034     if (name.Contains("friends")) continue;
1035
1036     // Get the path
1037     TString data(Form("%s/%s", file->GetTitle(), name.Data()));
1038
1039     TFile* test = TFile::Open(data.Data(), "READ");
1040     if (!test || test->IsZombie()) {
1041       ::Warning("ScanDirectory", "Failed to open file %s", data.Data());
1042       continue;
1043     }
1044     test->Close();
1045     chain->Add(data);
1046   }
1047 }
1048
1049 //____________________________________________________________________
1050 TChain*
1051 AliFMDInput::MakeChain(const char* what, const char* datadir, bool recursive)
1052 {
1053   TString w(what);
1054   w.ToUpper();
1055   const char* treeName = 0;
1056   const char* pattern  = 0;
1057   if      (w.Contains("ESD")) { treeName = "esdTree"; pattern = "AliESD"; }
1058   else if (w.Contains("MC"))  { treeName = "TE";      pattern = "galice"; }
1059   else {
1060     ::Error("MakeChain", "Unknown mode '%s' (not one of ESD, or MC)", what);
1061     return 0;
1062   }
1063
1064   // --- Our data chain ----------------------------------------------
1065   TChain* chain = new TChain(treeName);
1066
1067   // --- Get list of ESDs --------------------------------------------
1068   // Open source directory, and make sure we go back to were we were
1069   TString oldDir(gSystem->WorkingDirectory());
1070   TSystemDirectory d(datadir, datadir);
1071   ScanDirectory(&d, oldDir, chain, pattern, recursive);
1072
1073   return chain;
1074 }
1075
1076
1077 //____________________________________________________________________
1078 //
1079 // EOF
1080 //