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   delete ll;
 207 }
 208
 209
 210 //____________________________________________________________________
 211 AliFMDInput::ETrees
 212 AliFMDInput::ParseLoad(const char* what)
 213 {
 214   TString opt(what);
 215   opt.ToLower();
 216   const ETrees* ptype = fgkAllLoads;
 217   do {
 218     ETrees  type = *ptype;
 219     if (opt.Contains(TreeName(type,true), TString::kIgnoreCase))
 220       return type;
 221   } while (*ptype++ != kUser);
 222   return kUser;
 223 }
 224 //____________________________________________________________________
 225 const char*
 226 AliFMDInput::LoadedString(Bool_t dataOnly) const
 227 {
 228   static TString ret;
 229   if (!ret.IsNull()) return ret.Data();
 230
 231   const ETrees* ptype = fgkAllLoads;
 232   do {
 233     ETrees type = *ptype;
 234     if (dataOnly &&
 235         (type == kKinematics ||
 236          type == kHeader     ||
 237          type == kGeometry   ||
 238          type == kTrackRefs)) continue;
 239     if (!IsLoaded(*ptype)) continue;
 240
 241     if (!ret.IsNull()) ret.Append(",");
 242     ret.Append(TreeName(type));
 243   } while (*ptype++ != kUser);
 244   return ret.Data();
 245 }
 246
 247 //____________________________________________________________________
 248 const char*
 249 AliFMDInput::TreeName(ETrees tree, Bool_t shortest)
 250 {
 251   if (shortest) {
 252     switch (tree) {
 253     case kHits:         return "hit";
 254     case kKinematics:   return "kin";
 255     case kDigits:       return "dig";
 256     case kSDigits:      return "sdig";
 257     case kHeader:       return "hea";
 258     case kRecPoints:    return "recp";
 259     case kESD:          return "esd";
 260     case kRaw:          return "raw";
 261     case kGeometry:     return "geo";
 262     case kTrackRefs:    return "trackr";
 263     case kRawCalib:     return "rawc";
 264     case kUser:         return "user";
 265     }
 266     return 0;
 267   }
 268   switch (tree) {
 269   case kHits:           return "Hits";
 270   case kKinematics:     return "Kinematics";
 271   case kDigits:         return "Digits";
 272   case kSDigits:        return "SDigits";
 273   case kHeader:         return "Header";
 274   case kRecPoints:      return "RecPoints";
 275   case kESD:            return "ESD";
 276   case kRaw:            return "Raw";
 277   case kGeometry:       return "Geometry";
 278   case kTrackRefs:      return "TrackRefs";
 279   case kRawCalib:       return "RawCalib";
 280   case kUser:           return "User";
 281   }
 282   return 0;
 283 }
 284
 285 //____________________________________________________________________
 286 Int_t
 287 AliFMDInput::NEvents() const
 288 {
 289   // Get number of events
 290   if (IsLoaded(kRaw) ||
 291       IsLoaded(kRawCalib)) return fReader->GetNumberOfEvents();
 292   if (fChainE) return fChainE->GetEntriesFast();
 293   if (fTreeE) return fTreeE->GetEntries();
 294   return -1;
 295 }
 296
 297 //____________________________________________________________________
 298 Bool_t
 299 AliFMDInput::Init()
 300 {
 301   // Initialize the object.  Get the needed loaders, and such.
 302
 303   // Check if we have been initialized
 304   if (fIsInit) {
 305     AliWarning("Already initialized");
 306     return fIsInit;
 307   }
 308   TString what;
 309   const ETrees* ptype = fgkAllLoads;
 310   do {
 311     ETrees type = *ptype;
 312     what.Append(Form("\n\t%-20s: %s", TreeName(type),
 313                      IsLoaded(type) ? "yes" : "no"));
 314   } while (*ptype++ != kUser);
 315
 316   Info("Init","Initialising w/mask 0x%04x%s", fTreeMask, what.Data());
 317   // Get the run
 318   if (IsLoaded(kDigits)     ||
 319       IsLoaded(kSDigits)    ||
 320       IsLoaded(kKinematics) ||
 321       IsLoaded(kTrackRefs)  ||
 322       IsLoaded(kHeader)) {
 323     if (!gSystem->FindFile(".:/", fGAliceFile)) {
 324       AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
 325     }
 326     else {
 327       fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
 328       if (!fLoader) {
 329         AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
 330         return kFALSE;
 331       }
 332       AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
 333
 334       if  (fLoader->LoadgAlice()) return kFALSE;
 335
 336       fRun = fLoader->GetAliRun();
 337
 338       // Get the FMD
 339       fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
 340       if (!fFMD) {
 341         AliError("Failed to get detector FMD from loader");
 342         return kFALSE;
 343       }
 344
 345       // Get the FMD loader
 346       fFMDLoader = fLoader->GetLoader("FMDLoader");
 347       if (!fFMDLoader) {
 348         AliError("Failed to get detector FMD loader from loader");
 349         return kFALSE;
 350       }
 351       if (fLoader->LoadHeader()) {
 352         AliError("Failed to get event header information from loader");
 353         return kFALSE;
 354       }
 355       fTreeE = fLoader->TreeE();
 356     }
 357   }
 358
 359   // Optionally, get the ESD files
 360   if (IsLoaded(kESD)) {
 361     fChainE = MakeChain("ESD", fInputDir, true);
 362     fESDEvent = new AliESDEvent();
 363     fESDEvent->ReadFromTree(fChainE);
 364     //    fChainE->SetBranchAddress("ESD", &fMainESD);
 365
 366   }
 367
 368   if (IsLoaded(kRaw) ||
 369       IsLoaded(kRawCalib)) {
 370     AliInfo("Getting FMD raw data digits");
 371     fArrayA = new TClonesArray("AliFMDDigit");
 372 #if 0
 373     if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
 374       fReader = new AliRawReaderRoot(fRawFile.Data());
 375     else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
 376       fReader = new AliRawReaderDate(fRawFile.Data());
 377     else
 378       fReader = new AliRawReaderFile(-1);
 379 #else
 380     if(!fRawFile.IsNull())
 381       fReader = AliRawReader::Create(fRawFile.Data());
 382     else
 383       fReader = new AliRawReaderFile(-1);
 384 #endif
 385     fFMDReader = new AliFMDRawReader(fReader, 0);
 386   }
 387
 388   // Optionally, get the geometry
 389   if (IsLoaded(kGeometry)) {
 390     TString fname;
 391     if (fRun) {
 392       fname = gSystem->DirName(fGAliceFile);
 393       fname.Append("/geometry.root");
 394     }
 395     if (!gSystem->AccessPathName(fname.Data()))
 396       fname = "";
 397     AliCDBManager* cdb   = AliCDBManager::Instance();
 398     if (!cdb->IsDefaultStorageSet()) {
 399       cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 400       cdb->SetRun(0);
 401     }
 402
 403     AliGeomManager::LoadGeometry(fname.IsNull() ? 0 : fname.Data());
 404
 405     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
 406     if (align) {
 407       AliInfo("Got alignment data from CDB");
 408       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
 409       if (!array) {
 410         AliWarning("Invalid align data from CDB");
 411       }
 412       else {
 413         Int_t nAlign = array->GetEntries();
 414         for (Int_t i = 0; i < nAlign; i++) {
 415           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
 416           if (!a->ApplyToGeometry()) {
 417             AliWarning(Form("Failed to apply alignment to %s",
 418                             a->GetSymName()));
 419           }
 420         }
 421       }
 422     }
 423     AliFMDGeometry* geom = AliFMDGeometry::Instance();
 424     geom->Init();
 425     geom->InitTransformations();
 426   }
 427
 428   fEventCount = 0;
 429   fIsInit = kTRUE;
 430   return fIsInit;
 431 }
 432
 433 //____________________________________________________________________
 434 Bool_t
 435 AliFMDInput::Begin(Int_t event)
 436 {
 437   // Called at the begining of each event.  Per default, it gets the
 438   // data trees and gets pointers to the output arrays.   Users can
 439   // overload this, but should call this member function in the
 440   // overloaded member function of the derived class.
 441
 442   // Check if we have been initialized
 443   if (!fIsInit) {
 444     AliError("Not initialized");
 445     return fIsInit;
 446   }
 447
 448   // Get the event
 449   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
 450
 451   // Possibly load global kinematics information
 452   if (IsLoaded(kKinematics)) {
 453     // AliInfo("Getting kinematics");
 454     if (fLoader->LoadKinematics("READ")) return kFALSE;
 455     fStack = fLoader->Stack();
 456   }
 457
 458   // Possibly load FMD Hit information
 459   if (IsLoaded(kHits)) {
 460     // AliInfo("Getting FMD hits");
 461     if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
 462     fTreeH = fFMDLoader->TreeH();
 463     if (!fArrayH) fArrayH = fFMD->Hits();
 464   }
 465
 466   // Possibly load FMD TrackReference information
 467   if (IsLoaded(kTrackRefs)) {
 468     // AliInfo("Getting FMD hits");
 469     if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
 470     fTreeTR = fLoader->TreeTR();
 471     if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
 472     fTreeTR->SetBranchAddress("TrackReferences",  &fArrayTR);
 473   }
 474
 475   // Possibly load heaedr information
 476   if (IsLoaded(kHeader)) {
 477     // AliInfo("Getting FMD hits");
 478     if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
 479     fHeader = fLoader->GetHeader();
 480   }
 481
 482   // Possibly load FMD Digit information
 483   if (IsLoaded(kDigits)) {
 484     // AliInfo("Getting FMD digits");
 485     if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
 486     fTreeD = fFMDLoader->TreeD();
 487     if (fTreeD) {
 488       if (!fArrayD) fArrayD = fFMD->Digits();
 489     }
 490     else {
 491       fArrayD = 0;
 492       AliWarning(Form("Failed to load FMD Digits"));
 493     }
 494   }
 495
 496   // Possibly load FMD Sdigit information
 497   if (IsLoaded(kSDigits)) {
 498     // AliInfo("Getting FMD summable digits");
 499     if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) {
 500       AliWarning("Failed to load SDigits!");
 501       return kFALSE;
 502     }
 503     fTreeS = fFMDLoader->TreeS();
 504     if (!fArrayS) fArrayS = fFMD->SDigits();
 505   }
 506
 507   // Possibly load FMD RecPoints information
 508   if (IsLoaded(kRecPoints)) {
 509     // AliInfo("Getting FMD reconstructed points");
 510     if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
 511     fTreeR = fFMDLoader->TreeR();
 512     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
 513     fTreeR->SetBranchAddress("FMD",  &fArrayR);
 514   }
 515
 516   // Possibly load FMD ESD information
 517   if (IsLoaded(kESD)) {
 518     // AliInfo("Getting FMD event summary data");
 519     Int_t read = fChainE->GetEntry(event);
 520     if (read <= 0) return kFALSE;
 521     fESD = fESDEvent->GetFMDData();
 522     if (!fESD) return kFALSE;
 523   }
 524
 525   // Possibly load FMD Digit information
 526   if (IsLoaded(kRaw) || IsLoaded(kRawCalib)) {
 527     Bool_t mon = fRawFile.Contains("mem://");
 528     // AliInfo("Getting FMD raw data digits");
 529     if (mon) std::cout << "Waiting for event ..." << std::flush;
 530     do {
 531       if (!fReader->NextEvent()) {
 532         if (mon) {
 533           gSystem->Sleep(3);
 534           continue;
 535         }
 536         return kFALSE;
 537       }
 538       UInt_t eventType = fReader->GetType();
 539       if(eventType == AliRawEventHeaderBase::kPhysicsEvent ||
 540          eventType == AliRawEventHeaderBase::kCalibrationEvent)
 541         break;
 542     } while (true);
 543     if (mon) std::cout << "got it" << std::endl;
 544     // AliFMDRawReader r(fReader, 0);
 545     fArrayA->Clear();
 546     fFMDReader->ReadAdcs(fArrayA);
 547     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
 548   }
 549   fEventCount++;
 550   return kTRUE;
 551 }
 552
 553
 554 //____________________________________________________________________
 555 Bool_t
 556 AliFMDInput::Event()
 557 {
 558   // Process one event.  The default implementation one or more of
 559   //
 560   //   -  ProcessHits       if the hits are loaded.
 561   //   -  ProcessDigits     if the digits are loaded.
 562   //   -  ProcessSDigits    if the sumbable digits are loaded.
 563   //   -  ProcessRecPoints  if the reconstructed points are loaded.
 564   //   -  ProcessESD        if the event summary data is loaded
 565   //
 566   if (IsLoaded(kHits))     if (!ProcessHits())      return kFALSE;
 567   if (IsLoaded(kTrackRefs))if (!ProcessTrackRefs()) return kFALSE;
 568   if (IsLoaded(kKinematics) &&
 569       IsLoaded(kHits))     if (!ProcessTracks())    return kFALSE;
 570   if (IsLoaded(kKinematics))if (!ProcessStack())     return kFALSE;
 571   if (IsLoaded(kSDigits))  if (!ProcessSDigits())   return kFALSE;
 572   if (IsLoaded(kDigits))   if (!ProcessDigits())    return kFALSE;
 573   if (IsLoaded(kRaw))      if (!ProcessRawDigits()) return kFALSE;
 574   if (IsLoaded(kRawCalib)) if (!ProcessRawCalibDigits())return kFALSE;
 575   if (IsLoaded(kRecPoints))if (!ProcessRecPoints()) return kFALSE;
 576   if (IsLoaded(kESD))      if (!ProcessESDs())      return kFALSE;
 577   if (IsLoaded(kUser))     if (!ProcessUsers())     return kFALSE;
 578
 579   return kTRUE;
 580 }
 581
 582 //____________________________________________________________________
 583 Bool_t
 584 AliFMDInput::ProcessHits()
 585 {
 586   // Read the hit tree, and pass each hit to the member function
 587   // ProcessHit.
 588   if (!fTreeH) {
 589     AliError("No hit tree defined");
 590     return kFALSE;
 591   }
 592   if (!fArrayH) {
 593     AliError("No hit array defined");
 594     return kFALSE;
 595   }
 596
 597   Int_t nTracks = fTreeH->GetEntries();
 598   for (Int_t i = 0; i < nTracks; i++) {
 599     Int_t hitRead  = fTreeH->GetEntry(i);
 600     if (hitRead <= 0) continue;
 601
 602     Int_t nHit = fArrayH->GetEntries();
 603     if (nHit <= 0) continue;
 604
 605     for (Int_t j = 0; j < nHit; j++) {
 606       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 607       if (!hit) continue;
 608
 609       TParticle* track = 0;
 610       if (IsLoaded(kKinematics) && fStack) {
 611         Int_t trackno = hit->Track();
 612         track = fStack->Particle(trackno);
 613       }
 614       if (!ProcessHit(hit, track)) return kFALSE;
 615     }
 616   }
 617   return kTRUE;
 618 }
 619 //____________________________________________________________________
 620 Bool_t
 621 AliFMDInput::ProcessTrackRefs()
 622 {
 623   // Read the reconstrcted points tree, and pass each reconstruction
 624   // object (AliFMDRecPoint) to either ProcessRecPoint.
 625   if (!fTreeTR) {
 626     AliError("No track reference tree defined");
 627     return kFALSE;
 628   }
 629   if (!fArrayTR) {
 630     AliError("No track reference array defined");
 631     return kFALSE;
 632   }
 633
 634   Int_t nEv = fTreeTR->GetEntries();
 635   for (Int_t i = 0; i < nEv; i++) {
 636     Int_t trRead  = fTreeTR->GetEntry(i);
 637     if (trRead <= 0) continue;
 638     Int_t nTrackRefs = fArrayTR->GetEntries();
 639     for (Int_t j = 0; j < nTrackRefs; j++) {
 640       AliTrackReference* trackRef =
 641         static_cast<AliTrackReference*>(fArrayTR->At(j));
 642       if (!trackRef) continue;
 643       // if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
 644       TParticle* track = 0;
 645       if (IsLoaded(kKinematics) && fStack) {
 646         Int_t trackno = trackRef->GetTrack();
 647         track = fStack->Particle(trackno);
 648       }
 649       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
 650     }
 651   }
 652   return kTRUE;
 653 }
 654 //____________________________________________________________________
 655 Bool_t
 656 AliFMDInput::ProcessTracks()
 657 {
 658   // Read the hit tree, and pass each hit to the member function
 659   // ProcessTrack.
 660   if (!fStack) {
 661     AliError("No track tree defined");
 662     return kFALSE;
 663   }
 664   if (!fTreeH) {
 665     AliError("No hit tree defined");
 666     return kFALSE;
 667   }
 668   if (!fArrayH) {
 669     AliError("No hit array defined");
 670     return kFALSE;
 671   }
 672
 673   // Int_t nTracks = fStack->GetNtrack();
 674   Int_t nTracks = fTreeH->GetEntries();
 675   for (Int_t i = 0; i < nTracks; i++) {
 676     Int_t      trackno = nTracks - i - 1;
 677     TParticle* track   = fStack->Particle(trackno);
 678     if (!track) continue;
 679
 680     // Get the hits for this track.
 681     Int_t hitRead  = fTreeH->GetEntry(i);
 682     Int_t nHit     = fArrayH->GetEntries();
 683     if (nHit == 0 || hitRead <= 0) {
 684       // Let user code see the track, even if there's no hits.
 685       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
 686       continue;
 687     }
 688
 689     // Loop over the hits corresponding to this track.
 690     for (Int_t j = 0; j < nHit; j++) {
 691       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 692       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
 693     }
 694   }
 695   return kTRUE;
 696 }
 697 //____________________________________________________________________
 698 Bool_t
 699 AliFMDInput::ProcessStack()
 700 {
 701   // Read the hit tree, and pass each hit to the member function
 702   // ProcessTrack.
 703   if (!fStack) {
 704     AliError("No track tree defined");
 705     return kFALSE;
 706   }
 707   Int_t nTracks = fStack->GetNtrack();
 708   for (Int_t i = 0; i < nTracks; i++) {
 709     Int_t      trackno = nTracks - i - 1;
 710     TParticle* track   = fStack->Particle(trackno);
 711     if (!track) continue;
 712
 713     if (!ProcessParticle(trackno, track)) return kFALSE;
 714   }
 715   return kTRUE;
 716 }
 717 //____________________________________________________________________
 718 Bool_t
 719 AliFMDInput::ProcessDigits()
 720 {
 721   // Read the digit tree, and pass each digit to the member function
 722   // ProcessDigit.
 723   if (!fTreeD) {
 724     AliError("No digit tree defined");
 725     return kFALSE;
 726   }
 727   if (!fArrayD) {
 728     AliError("No digit array defined");
 729     return kFALSE;
 730   }
 731
 732   Int_t nEv = fTreeD->GetEntries();
 733   for (Int_t i = 0; i < nEv; i++) {
 734     Int_t digitRead  = fTreeD->GetEntry(i);
 735     if (digitRead <= 0) continue;
 736     Int_t nDigit = fArrayD->GetEntries();
 737     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
 738     if (nDigit <= 0) continue;
 739     for (Int_t j = 0; j < nDigit; j++) {
 740       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 741       if (!digit) continue;
 742       if (!ProcessDigit(digit)) return kFALSE;
 743     }
 744   }
 745   return kTRUE;
 746 }
 747
 748 //____________________________________________________________________
 749 Bool_t
 750 AliFMDInput::ProcessSDigits()
 751 {
 752   // Read the summable digit tree, and pass each sumable digit to the
 753   // member function ProcessSdigit.
 754   if (!fTreeS) {
 755     AliWarning("No sdigit tree defined");
 756     return kTRUE; // Empty SDigits is fine
 757   }
 758   if (!fArrayS) {
 759     AliWarning("No sdigit array defined");
 760     return kTRUE; // Empty SDigits is fine
 761   }
 762
 763   Int_t nEv = fTreeS->GetEntries();
 764   for (Int_t i = 0; i < nEv; i++) {
 765     Int_t sdigitRead  = fTreeS->GetEntry(i);
 766     if (sdigitRead <= 0) {
 767       AliInfo(Form("Read nothing from tree"));
 768       continue;
 769     }
 770     Int_t nSdigit = fArrayS->GetEntriesFast();
 771     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
 772     AliInfo(Form("Got %5d digits for this event", nSdigit));
 773     if (nSdigit <= 0) continue;
 774     for (Int_t j = 0; j < nSdigit; j++) {
 775       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 776       if (!sdigit) continue;
 777       if (!ProcessSDigit(sdigit)) return kFALSE;
 778     }
 779   }
 780   return kTRUE;
 781 }
 782
 783 //____________________________________________________________________
 784 Bool_t
 785 AliFMDInput::ProcessRawDigits()
 786 {
 787   // Read the digit tree, and pass each digit to the member function
 788   // ProcessDigit.
 789   if (!fArrayA) {
 790     AliError("No raw digit array defined");
 791     return kFALSE;
 792   }
 793
 794   Int_t nDigit = fArrayA->GetEntries();
 795   if (nDigit <= 0) return kTRUE;
 796   for (Int_t j = 0; j < nDigit; j++) {
 797     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 798     if (!digit) continue;
 799     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 800       digit->Print();
 801     if (!ProcessRawDigit(digit)) return kFALSE;
 802   }
 803   return kTRUE;
 804 }
 805
 806 //____________________________________________________________________
 807 Bool_t
 808 AliFMDInput::ProcessRawCalibDigits()
 809 {
 810   // Read the digit tree, and pass each digit to the member function
 811   // ProcessDigit.
 812   if (!fArrayA) {
 813     AliError("No raw digit array defined");
 814     return kFALSE;
 815   }
 816
 817   Int_t nDigit = fArrayA->GetEntries();
 818   if (nDigit <= 0) return kTRUE;
 819   for (Int_t j = 0; j < nDigit; j++) {
 820     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 821     if (!digit) continue;
 822     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 823       digit->Print();
 824     if (!ProcessRawCalibDigit(digit)) return kFALSE;
 825   }
 826   return kTRUE;
 827 }
 828
 829 //____________________________________________________________________
 830 Bool_t
 831 AliFMDInput::ProcessRecPoints()
 832 {
 833   // Read the reconstrcted points tree, and pass each reconstruction
 834   // object (AliFMDRecPoint) to either ProcessRecPoint.
 835   if (!fTreeR) {
 836     AliError("No recpoint tree defined");
 837     return kFALSE;
 838   }
 839   if (!fArrayR) {
 840     AliError("No recpoints array defined");
 841     return kFALSE;
 842   }
 843
 844   Int_t nEv = fTreeR->GetEntries();
 845   for (Int_t i = 0; i < nEv; i++) {
 846     Int_t recRead  = fTreeR->GetEntry(i);
 847     if (recRead <= 0) continue;
 848     Int_t nRecPoint = fArrayR->GetEntries();
 849     for (Int_t j = 0; j < nRecPoint; j++) {
 850       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 851       if (!recPoint) continue;
 852       if (!ProcessRecPoint(recPoint)) return kFALSE;
 853     }
 854   }
 855   return kTRUE;
 856 }
 857
 858 //____________________________________________________________________
 859 Bool_t
 860 AliFMDInput::ProcessESDs()
 861 {
 862   // Process event summary data
 863   if (!fESD) return kFALSE;
 864   for (UShort_t det = 1; det <= 3; det++) {
 865     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 866     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 867       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 868       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 869       for (UShort_t sec = 0; sec < nsec; sec++) {
 870         for (UShort_t str = 0; str < nstr; str++) {
 871           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 872           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 873           if (!fESD->IsAngleCorrected())
 874             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 875           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 876         }
 877       }
 878     }
 879   }
 880   return kTRUE;
 881 }
 882
 883 //____________________________________________________________________
 884 Bool_t
 885 AliFMDInput::ProcessUsers()
 886 {
 887   // Process event summary data
 888   for (UShort_t det = 1; det <= 3; det++) {
 889     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 890     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 891       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 892       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 893       for (UShort_t sec = 0; sec < nsec; sec++) {
 894         for (UShort_t str = 0; str < nstr; str++) {
 895           Float_t v  = GetSignal(det,*rng,sec,str);
 896           if (!ProcessUser(det, *rng, sec, str, v)) continue;
 897         }
 898       }
 899     }
 900   }
 901   return kTRUE;
 902 }
 903
 904 //____________________________________________________________________
 905 Bool_t
 906 AliFMDInput::End()
 907 {
 908   // Called at the end of each event.  Per default, it unloads the
 909   // data trees and resets the pointers to the output arrays.   Users
 910   // can overload this, but should call this member function in the
 911   // overloaded member function of the derived class.
 912
 913   // Check if we have been initialized
 914   if (!fIsInit) {
 915     AliError("Not initialized");
 916     return fIsInit;
 917   }
 918   // Possibly unload global kinematics information
 919   if (IsLoaded(kKinematics)) {
 920     fLoader->UnloadKinematics();
 921     // fTreeK = 0;
 922     fStack = 0;
 923   }
 924   // Possibly unload FMD Hit information
 925   if (IsLoaded(kHits)) {
 926     fFMDLoader->UnloadHits();
 927     fTreeH = 0;
 928   }
 929   // Possibly unload FMD Digit information
 930   if (IsLoaded(kDigits)) {
 931     fFMDLoader->UnloadDigits();
 932     fTreeD = 0;
 933   }
 934   // Possibly unload FMD Sdigit information
 935   if (IsLoaded(kSDigits)) {
 936     fFMDLoader->UnloadSDigits();
 937     fTreeS = 0;
 938   }
 939   // Possibly unload FMD RecPoints information
 940   if (IsLoaded(kRecPoints)) {
 941     fFMDLoader->UnloadRecPoints();
 942     fTreeR = 0;
 943   }
 944   // AliInfo("Now out event");
 945   return kTRUE;
 946 }
 947
 948 //____________________________________________________________________
 949 Bool_t
 950 AliFMDInput::Run(UInt_t maxEvents)
 951 {
 952   // Run over all events and files references in galice.root
 953
 954   Bool_t retval;
 955   if (!(retval = Init())) return retval;
 956
 957   fNEvents = NEvents();
 958   if (fNEvents < 0)       fNEvents = maxEvents;
 959   else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
 960
 961   Int_t event = 0;
 962   for (; fNEvents < 0 || event < fNEvents; event++) {
 963     printf("\rEvent %8d/%8d ...", event, fNEvents);
 964     if (!(retval = Begin(event))) break;
 965     if (!(retval = Event())) break;
 966     if (!(retval = End())) break;
 967   }
 968   printf("Looped over %8d events\n", event+1);
 969   if (!retval) return retval;
 970   retval = Finish();
 971   return retval;
 972 }
 973
 974 //__________________________________________________________________
 975 TArrayF
 976 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
 977 {
 978   // Service function to define a logarithmic axis.
 979   // Parameters:
 980   //   n    Number of bins
 981   //   min  Minimum of axis
 982   //   max  Maximum of axis
 983   TArrayF bins(n+1);
 984   bins[0]      = min;
 985   if (n <= 20) {
 986     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
 987     return bins;
 988   }
 989   Float_t dp   = n / TMath::Log10(max / min);
 990   Float_t pmin = TMath::Log10(min);
 991   for (Int_t i = 1; i < n+1; i++) {
 992     Float_t p = pmin + i / dp;
 993     bins[i]   = TMath::Power(10, p);
 994   }
 995   return bins;
 996 }
 997
 998 //____________________________________________________________________
 999 void
1000 AliFMDInput::ScanDirectory(TSystemDirectory* dir,
1001                            const TString& olddir,
1002                            TChain* chain,
1003                            const char* pattern, bool recursive)
1004 {
1005   // Get list of files, and go back to old working directory
1006   TString oldDir(gSystem->WorkingDirectory());
1007   TList* files = dir->GetListOfFiles();
1008   gSystem->ChangeDirectory(oldDir);
1009
1010   // Sort list of files and check if we should add it
1011   if (!files) return;
1012   files->Sort();
1013   TIter next(files);
1014   TSystemFile* file = 0;
1015   while ((file = static_cast<TSystemFile*>(next()))) {
1016     TString name(file->GetName());
1017
1018     // Ignore special links
1019     if (name == "." || name == "..") continue;
1020
1021     // Check if this is a directory
1022     if (file->IsDirectory()) {
1023       if (recursive)
1024         ScanDirectory(static_cast<TSystemDirectory*>(file),
1025                       olddir, chain,
1026                       pattern,recursive);
1027       continue;
1028     }
1029
1030     // If this is not a root file, ignore
1031     if (!name.EndsWith(".root")) continue;
1032
1033     // If this file does not contain the pattern, ignore
1034     if (!name.Contains(pattern)) continue;
1035     if (name.Contains("friends")) continue;
1036
1037     // Get the path
1038     TString data(Form("%s/%s", file->GetTitle(), name.Data()));
1039
1040     TFile* test = TFile::Open(data.Data(), "READ");
1041     if (!test || test->IsZombie()) {
1042       ::Warning("ScanDirectory", "Failed to open file %s", data.Data());
1043       continue;
1044     }
1045     test->Close();
1046     chain->Add(data);
1047   }
1048 }
1049
1050 //____________________________________________________________________
1051 TChain*
1052 AliFMDInput::MakeChain(const char* what, const char* datadir, bool recursive)
1053 {
1054   TString w(what);
1055   w.ToUpper();
1056   const char* treeName = 0;
1057   const char* pattern  = 0;
1058   if      (w.Contains("ESD")) { treeName = "esdTree"; pattern = "AliESD"; }
1059   else if (w.Contains("MC"))  { treeName = "TE";      pattern = "galice"; }
1060   else {
1061     ::Error("MakeChain", "Unknown mode '%s' (not one of ESD, or MC)", what);
1062     return 0;
1063   }
1064
1065   // --- Our data chain ----------------------------------------------
1066   TChain* chain = new TChain(treeName);
1067
1068   // --- Get list of ESDs --------------------------------------------
1069   // Open source directory, and make sure we go back to were we were
1070   TString oldDir(gSystem->WorkingDirectory());
1071   TSystemDirectory d(datadir, datadir);
1072   ScanDirectory(&d, oldDir, chain, pattern, recursive);
1073
1074   return chain;
1075 }
1076
1077
1078 //____________________________________________________________________
1079 //
1080 // EOF
1081 //