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 "AliFMD.h"             // ALIFMD_H
  41 #include "AliFMDHit.h"          // ALIFMDHIT_H
  42 #include "AliFMDDigit.h"        // ALIFMDDigit_H
  43 #include "AliFMDSDigit.h"       // ALIFMDDigit_H
  44 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
  45 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
  46 #include <AliESD.h>
  47 #include <AliESDFMD.h>
  48 #include <AliESDEvent.h>
  49 #include <AliCDBManager.h>
  50 #include <AliCDBEntry.h>
  51 #include <AliAlignObjParams.h>
  52 #include <AliTrackReference.h>
  53 #include <TTree.h>              // ROOT_TTree
  54 #include <TChain.h>             // ROOT_TChain
  55 #include <TParticle.h>          // ROOT_TParticle
  56 #include <TString.h>            // ROOT_TString
  57 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
  58 #include <TMath.h>              // ROOT_TMath
  59 #include <TGeoManager.h>        // ROOT_TGeoManager
  60 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
  61 #include <Riostream.h>          // ROOT_Riostream
  62 #include <TFile.h>              // ROOT_TFile
  63 #include <TStreamerInfo.h>
  64 #include <TArrayF.h>
  65
  66 //____________________________________________________________________
  67 ClassImp(AliFMDInput)
  68 #if 0
  69   ; // This is here to keep Emacs for indenting the next line
  70 #endif
  71
  72
  73 //____________________________________________________________________
  74 AliFMDInput::AliFMDInput()
  75   : TNamed("AliFMDInput", "Input handler for various FMD data"),
  76     fGAliceFile(""),
  77     fLoader(0),
  78     fRun(0),
  79     fStack(0),
  80     fFMDLoader(0),
  81     fReader(0),
  82     fFMDReader(0),
  83     fFMD(0),
  84     fESD(0),
  85     fESDEvent(0),
  86     fTreeE(0),
  87     fTreeH(0),
  88     fTreeTR(0),
  89     fTreeD(0),
  90     fTreeS(0),
  91     fTreeR(0),
  92     fTreeA(0),
  93     fChainE(0),
  94     fArrayE(0),
  95     fArrayH(0),
  96     fArrayTR(0),
  97     fArrayD(0),
  98     fArrayS(0),
  99     fArrayR(0),
 100     fArrayA(0),
 101     fHeader(0),
 102     fGeoManager(0),
 103     fTreeMask(0),
 104     fRawFile(""),
 105     fIsInit(kFALSE),
 106     fEventCount(0)
 107 {
 108
 109   // Constructor of an FMD input object.  Specify what data to read in
 110   // using the AddLoad member function.   Sub-classes should at a
 111   // minimum overload the member function Event.   A full job can be
 112   // executed using the member function Run.
 113 }
 114
 115
 116
 117 //____________________________________________________________________
 118 AliFMDInput::AliFMDInput(const char* gAliceFile)
 119   : TNamed("AliFMDInput", "Input handler for various FMD data"),
 120     fGAliceFile(gAliceFile),
 121     fLoader(0),
 122     fRun(0),
 123     fStack(0),
 124     fFMDLoader(0),
 125     fReader(0),
 126     fFMDReader(0),
 127     fFMD(0),
 128     fESD(0),
 129     fESDEvent(0),
 130     fTreeE(0),
 131     fTreeH(0),
 132     fTreeTR(0),
 133     fTreeD(0),
 134     fTreeS(0),
 135     fTreeR(0),
 136     fTreeA(0),
 137     fChainE(0),
 138     fArrayE(0),
 139     fArrayH(0),
 140     fArrayTR(0),
 141     fArrayD(0),
 142     fArrayS(0),
 143     fArrayR(0),
 144     fArrayA(0),
 145     fHeader(0),
 146     fGeoManager(0),
 147     fTreeMask(0),
 148     fRawFile(""),
 149     fIsInit(kFALSE),
 150     fEventCount(0)
 151 {
 152
 153   // Constructor of an FMD input object.  Specify what data to read in
 154   // using the AddLoad member function.   Sub-classes should at a
 155   // minimum overload the member function Event.   A full job can be
 156   // executed using the member function Run.
 157 }
 158
 159 //____________________________________________________________________
 160 Int_t
 161 AliFMDInput::NEvents() const
 162 {
 163   // Get number of events
 164   if (TESTBIT(fTreeMask, kRaw)) return fReader->GetNumberOfEvents();
 165   if (fChainE) return fChainE->GetEntriesFast();
 166   if (fTreeE) return fTreeE->GetEntries();
 167   return -1;
 168 }
 169
 170 //____________________________________________________________________
 171 Bool_t
 172 AliFMDInput::Init()
 173 {
 174   // Initialize the object.  Get the needed loaders, and such.
 175
 176   // Check if we have been initialized
 177   if (fIsInit) {
 178     AliWarning("Already initialized");
 179     return fIsInit;
 180   }
 181   Info("Init","Initialising w/mask 0x%04x\n"
 182        "\tHits:          %d\n"
 183        "\tKinematics:    %d\n"
 184        "\tDigits:        %d\n"
 185        "\tSDigits:       %d\n"
 186        "\tHeader:        %d\n"
 187        "\tRecPoints:     %d\n"
 188        "\tESD:           %d\n"
 189        "\tRaw:           %d\n"
 190        "\tGeometry:      %d\n"
 191        "\tTracks:        %d\n"
 192        "\tTracksRefs:    %d",
 193        fTreeMask,
 194        TESTBIT(fTreeMask, kHits),
 195        TESTBIT(fTreeMask, kKinematics),
 196        TESTBIT(fTreeMask, kDigits),
 197        TESTBIT(fTreeMask, kSDigits),
 198        TESTBIT(fTreeMask, kHeader),
 199        TESTBIT(fTreeMask, kRecPoints),
 200        TESTBIT(fTreeMask, kESD),
 201        TESTBIT(fTreeMask, kRaw),
 202        TESTBIT(fTreeMask, kGeometry),
 203        TESTBIT(fTreeMask, kTracks),
 204        TESTBIT(fTreeMask, kTrackRefs));
 205   // Get the run
 206   if (TESTBIT(fTreeMask, kDigits)     ||
 207       TESTBIT(fTreeMask, kSDigits)    ||
 208       TESTBIT(fTreeMask, kKinematics) ||
 209       TESTBIT(fTreeMask, kTrackRefs)  ||
 210       TESTBIT(fTreeMask, kTracks)     ||
 211       TESTBIT(fTreeMask, kHeader)) {
 212     if (!gSystem->FindFile(".:/", fGAliceFile)) {
 213       AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
 214     }
 215     else {
 216       fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
 217       if (!fLoader) {
 218         AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
 219         return kFALSE;
 220       }
 221       AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
 222
 223       if  (fLoader->LoadgAlice()) return kFALSE;
 224
 225       fRun = fLoader->GetAliRun();
 226
 227       // Get the FMD
 228       fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
 229       if (!fFMD) {
 230         AliError("Failed to get detector FMD from loader");
 231         return kFALSE;
 232       }
 233
 234       // Get the FMD loader
 235       fFMDLoader = fLoader->GetLoader("FMDLoader");
 236       if (!fFMDLoader) {
 237         AliError("Failed to get detector FMD loader from loader");
 238         return kFALSE;
 239       }
 240       if (fLoader->LoadHeader()) {
 241         AliError("Failed to get event header information from loader");
 242         return kFALSE;
 243       }
 244       fTreeE = fLoader->TreeE();
 245     }
 246   }
 247
 248   // Optionally, get the ESD files
 249   if (TESTBIT(fTreeMask, kESD)) {
 250     fChainE = new TChain("esdTree");
 251     TSystemDirectory dir(".",".");
 252     TList*           files = dir.GetListOfFiles();
 253     TSystemFile*     file = 0;
 254     if (!files) {
 255       AliError("No files");
 256       return kFALSE;
 257     }
 258     files->Sort();
 259     TIter            next(files);
 260     while ((file = static_cast<TSystemFile*>(next()))) {
 261       TString fname(file->GetName());
 262       if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
 263     }
 264     fESDEvent = new AliESDEvent();
 265     fESDEvent->ReadFromTree(fChainE);
 266     //    fChainE->SetBranchAddress("ESD", &fMainESD);
 267
 268   }
 269
 270   if (TESTBIT(fTreeMask, kRaw)) {
 271     AliInfo("Getting FMD raw data digits");
 272     fArrayA = new TClonesArray("AliFMDDigit");
 273 #if 0
 274     if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
 275       fReader = new AliRawReaderRoot(fRawFile.Data());
 276     else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
 277       fReader = new AliRawReaderDate(fRawFile.Data());
 278     else
 279       fReader = new AliRawReaderFile(-1);
 280 #else
 281     if(!fRawFile.IsNull())
 282       fReader = AliRawReader::Create(fRawFile.Data());
 283     else
 284       fReader = new AliRawReaderFile(-1);
 285 #endif
 286     fFMDReader = new AliFMDRawReader(fReader, 0);
 287   }
 288
 289   // Optionally, get the geometry
 290   if (TESTBIT(fTreeMask, kGeometry)) {
 291     if (fRun) {
 292       TString fname(fRun->GetGeometryFileName());
 293       if (fname.IsNull()) {
 294         Warning("Init", "No file name for the geometry from AliRun");
 295         fname = gSystem->DirName(fGAliceFile);
 296         fname.Append("/geometry.root");
 297       }
 298       fGeoManager = TGeoManager::Import(fname.Data());
 299       if (!fGeoManager) {
 300         Fatal("Init", "No geometry manager found");
 301         return kFALSE;
 302       }
 303     }
 304     else {
 305       AliGeomManager::LoadGeometry();
 306     }
 307     AliCDBManager* cdb   = AliCDBManager::Instance();
 308     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
 309     if (align) {
 310       AliInfo("Got alignment data from CDB");
 311       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
 312       if (!array) {
 313         AliWarning("Invalid align data from CDB");
 314       }
 315       else {
 316         Int_t nAlign = array->GetEntries();
 317         for (Int_t i = 0; i < nAlign; i++) {
 318           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
 319           if (!a->ApplyToGeometry()) {
 320             AliWarning(Form("Failed to apply alignment to %s",
 321                             a->GetSymName()));
 322           }
 323         }
 324       }
 325     }
 326   }
 327
 328   fEventCount = 0;
 329   fIsInit = kTRUE;
 330   return fIsInit;
 331 }
 332
 333 //____________________________________________________________________
 334 Bool_t
 335 AliFMDInput::Begin(Int_t event)
 336 {
 337   // Called at the begining of each event.  Per default, it gets the
 338   // data trees and gets pointers to the output arrays.   Users can
 339   // overload this, but should call this member function in the
 340   // overloaded member function of the derived class.
 341
 342   // Check if we have been initialized
 343   if (!fIsInit) {
 344     AliError("Not initialized");
 345     return fIsInit;
 346   }
 347
 348   // Get the event
 349   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
 350   AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
 351
 352   // Possibly load global kinematics information
 353   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
 354     // AliInfo("Getting kinematics");
 355     if (fLoader->LoadKinematics("READ")) return kFALSE;
 356     fStack = fLoader->Stack();
 357   }
 358
 359   // Possibly load FMD Hit information
 360   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
 361     // AliInfo("Getting FMD hits");
 362     if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
 363     fTreeH = fFMDLoader->TreeH();
 364     if (!fArrayH) fArrayH = fFMD->Hits();
 365   }
 366
 367   // Possibly load FMD TrackReference information
 368   if (TESTBIT(fTreeMask, kTrackRefs) || TESTBIT(fTreeMask, kTracks)) {
 369     // AliInfo("Getting FMD hits");
 370     if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
 371     fTreeTR = fLoader->TreeTR();
 372     if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
 373     fTreeTR->SetBranchAddress("TrackReferences",  &fArrayTR);
 374   }
 375
 376   // Possibly load heaedr information
 377   if (TESTBIT(fTreeMask, kHeader)) {
 378     // AliInfo("Getting FMD hits");
 379     if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
 380     fHeader = fLoader->GetHeader();
 381   }
 382
 383   // Possibly load FMD Digit information
 384   if (TESTBIT(fTreeMask, kDigits)) {
 385     // AliInfo("Getting FMD digits");
 386     if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
 387     fTreeD = fFMDLoader->TreeD();
 388     if (fTreeD) {
 389       if (!fArrayD) fArrayD = fFMD->Digits();
 390     }
 391     else {
 392       fArrayD = 0;
 393       AliWarning(Form("Failed to load FMD Digits"));
 394     }
 395   }
 396
 397   // Possibly load FMD Sdigit information
 398   if (TESTBIT(fTreeMask, kSDigits)) {
 399     // AliInfo("Getting FMD summable digits");
 400     if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) {
 401       AliWarning("Failed to load SDigits!");
 402       return kFALSE;
 403     }
 404     fTreeS = fFMDLoader->TreeS();
 405     if (!fArrayS) fArrayS = fFMD->SDigits();
 406   }
 407
 408   // Possibly load FMD RecPoints information
 409   if (TESTBIT(fTreeMask, kRecPoints)) {
 410     // AliInfo("Getting FMD reconstructed points");
 411     if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
 412     fTreeR = fFMDLoader->TreeR();
 413     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
 414     fTreeR->SetBranchAddress("FMD",  &fArrayR);
 415   }
 416
 417   // Possibly load FMD ESD information
 418   if (TESTBIT(fTreeMask, kESD)) {
 419     // AliInfo("Getting FMD event summary data");
 420     Int_t read = fChainE->GetEntry(event);
 421     if (read <= 0) return kFALSE;
 422     fESD = fESDEvent->GetFMDData();
 423     if (!fESD) return kFALSE;
 424 #if 0
 425     TFile* f = fChainE->GetFile();
 426     if (f) {
 427       TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
 428       if (o) {
 429         TStreamerInfo* info = static_cast<TStreamerInfo*>(o);
 430         std::cout << "AliFMDMap class version read is "
 431                   <<  info->GetClassVersion() << std::endl;
 432       }
 433     }
 434     // fESD->CheckNeedUShort(fChainE->GetFile());
 435 #endif
 436   }
 437
 438   // Possibly load FMD Digit information
 439   if (TESTBIT(fTreeMask, kRaw)) {
 440     // AliInfo("Getting FMD raw data digits");
 441     if (!fReader->NextEvent()) return kFALSE;
 442     // AliFMDRawReader r(fReader, 0);
 443     fArrayA->Clear();
 444     fFMDReader->ReadAdcs(fArrayA);
 445     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
 446   }
 447   fEventCount++;
 448   return kTRUE;
 449 }
 450
 451
 452 //____________________________________________________________________
 453 Bool_t
 454 AliFMDInput::Event()
 455 {
 456   // Process one event.  The default implementation one or more of
 457   //
 458   //   -  ProcessHits       if the hits are loaded.
 459   //   -  ProcessDigits     if the digits are loaded.
 460   //   -  ProcessSDigits    if the sumbable digits are loaded.
 461   //   -  ProcessRecPoints  if the reconstructed points are loaded.
 462   //   -  ProcessESD        if the event summary data is loaded
 463   //
 464   if (TESTBIT(fTreeMask, kHits))
 465     if (!ProcessHits()) return kFALSE;
 466   if (TESTBIT(fTreeMask, kTrackRefs))
 467     if (!ProcessTrackRefs()) return kFALSE;
 468   if (TESTBIT(fTreeMask, kTracks))
 469     if (!ProcessTracks()) return kFALSE;
 470   if (TESTBIT(fTreeMask, kSDigits))
 471     if (!ProcessSDigits()) return kFALSE;
 472   if (TESTBIT(fTreeMask, kDigits))
 473     if (!ProcessDigits()) return kFALSE;
 474   if (TESTBIT(fTreeMask, kRaw))
 475     if (!ProcessRawDigits()) return kFALSE;
 476   if (TESTBIT(fTreeMask, kRecPoints))
 477     if (!ProcessRecPoints()) return kFALSE;
 478   if (TESTBIT(fTreeMask, kESD))
 479     if (!ProcessESDs()) return kFALSE;
 480
 481   return kTRUE;
 482 }
 483
 484 //____________________________________________________________________
 485 Bool_t
 486 AliFMDInput::ProcessHits()
 487 {
 488   // Read the hit tree, and pass each hit to the member function
 489   // ProcessHit.
 490   if (!fTreeH) {
 491     AliError("No hit tree defined");
 492     return kFALSE;
 493   }
 494   if (!fArrayH) {
 495     AliError("No hit array defined");
 496     return kFALSE;
 497   }
 498
 499   Int_t nTracks = fTreeH->GetEntries();
 500   for (Int_t i = 0; i < nTracks; i++) {
 501     Int_t hitRead  = fTreeH->GetEntry(i);
 502     if (hitRead <= 0) continue;
 503
 504     Int_t nHit = fArrayH->GetEntries();
 505     if (nHit <= 0) continue;
 506
 507     for (Int_t j = 0; j < nHit; j++) {
 508       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 509       if (!hit) continue;
 510
 511       TParticle* track = 0;
 512       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 513         Int_t trackno = hit->Track();
 514         track = fStack->Particle(trackno);
 515       }
 516       if (!ProcessHit(hit, track)) return kFALSE;
 517     }
 518   }
 519   return kTRUE;
 520 }
 521 //____________________________________________________________________
 522 Bool_t
 523 AliFMDInput::ProcessTrackRefs()
 524 {
 525   // Read the reconstrcted points tree, and pass each reconstruction
 526   // object (AliFMDRecPoint) to either ProcessRecPoint.
 527   if (!fTreeTR) {
 528     AliError("No track reference tree defined");
 529     return kFALSE;
 530   }
 531   if (!fArrayTR) {
 532     AliError("No track reference array defined");
 533     return kFALSE;
 534   }
 535
 536   Int_t nEv = fTreeTR->GetEntries();
 537   for (Int_t i = 0; i < nEv; i++) {
 538     Int_t trRead  = fTreeTR->GetEntry(i);
 539     if (trRead <= 0) continue;
 540     Int_t nTrackRefs = fArrayTR->GetEntries();
 541     for (Int_t j = 0; j < nTrackRefs; j++) {
 542       AliTrackReference* trackRef = static_cast<AliTrackReference*>(fArrayTR->At(j));
 543       if (!trackRef) continue;
 544       if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
 545       TParticle* track = 0;
 546       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 547         Int_t trackno = trackRef->GetTrack();
 548         track = fStack->Particle(trackno);
 549       }
 550       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
 551     }
 552   }
 553   return kTRUE;
 554 }
 555 //____________________________________________________________________
 556 Bool_t
 557 AliFMDInput::ProcessTracks()
 558 {
 559   // Read the hit tree, and pass each hit to the member function
 560   // ProcessTrack.
 561   if (!fStack) {
 562     AliError("No track tree defined");
 563     return kFALSE;
 564   }
 565   if (!fTreeH) {
 566     AliError("No hit tree defined");
 567     return kFALSE;
 568   }
 569   if (!fArrayH) {
 570     AliError("No hit array defined");
 571     return kFALSE;
 572   }
 573
 574   // Int_t nTracks = fStack->GetNtrack();
 575   Int_t nTracks = fTreeH->GetEntries();
 576   for (Int_t i = 0; i < nTracks; i++) {
 577     Int_t      trackno = nTracks - i - 1;
 578     TParticle* track   = fStack->Particle(trackno);
 579     if (!track) continue;
 580
 581     // Get the hits for this track.
 582     Int_t hitRead  = fTreeH->GetEntry(i);
 583     Int_t nHit     = fArrayH->GetEntries();
 584     if (nHit == 0 || hitRead <= 0) {
 585       // Let user code see the track, even if there's no hits.
 586       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
 587       continue;
 588     }
 589
 590     // Loop over the hits corresponding to this track.
 591     for (Int_t j = 0; j < nHit; j++) {
 592       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 593       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
 594     }
 595   }
 596   return kTRUE;
 597 }
 598
 599 //____________________________________________________________________
 600 Bool_t
 601 AliFMDInput::ProcessDigits()
 602 {
 603   // Read the digit tree, and pass each digit to the member function
 604   // ProcessDigit.
 605   if (!fTreeD) {
 606     AliError("No digit tree defined");
 607     return kFALSE;
 608   }
 609   if (!fArrayD) {
 610     AliError("No digit array defined");
 611     return kFALSE;
 612   }
 613
 614   Int_t nEv = fTreeD->GetEntries();
 615   for (Int_t i = 0; i < nEv; i++) {
 616     Int_t digitRead  = fTreeD->GetEntry(i);
 617     if (digitRead <= 0) continue;
 618     Int_t nDigit = fArrayD->GetEntries();
 619     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
 620     if (nDigit <= 0) continue;
 621     for (Int_t j = 0; j < nDigit; j++) {
 622       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 623       if (!digit) continue;
 624       if (!ProcessDigit(digit)) return kFALSE;
 625     }
 626   }
 627   return kTRUE;
 628 }
 629
 630 //____________________________________________________________________
 631 Bool_t
 632 AliFMDInput::ProcessSDigits()
 633 {
 634   // Read the summable digit tree, and pass each sumable digit to the
 635   // member function ProcessSdigit.
 636   if (!fTreeS) {
 637     AliWarning("No sdigit tree defined");
 638     return kTRUE; // Empty SDigits is fine
 639   }
 640   if (!fArrayS) {
 641     AliWarning("No sdigit array defined");
 642     return kTRUE; // Empty SDigits is fine
 643   }
 644
 645   Int_t nEv = fTreeS->GetEntries();
 646   for (Int_t i = 0; i < nEv; i++) {
 647     Int_t sdigitRead  = fTreeS->GetEntry(i);
 648     if (sdigitRead <= 0) {
 649       AliInfo(Form("Read nothing from tree"));
 650       continue;
 651     }
 652     Int_t nSdigit = fArrayS->GetEntriesFast();
 653     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
 654     AliInfo(Form("Got %5d digits for this event", nSdigit));
 655     if (nSdigit <= 0) continue;
 656     for (Int_t j = 0; j < nSdigit; j++) {
 657       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 658       if (!sdigit) continue;
 659       if (!ProcessSDigit(sdigit)) return kFALSE;
 660     }
 661   }
 662   return kTRUE;
 663 }
 664
 665 //____________________________________________________________________
 666 Bool_t
 667 AliFMDInput::ProcessRawDigits()
 668 {
 669   // Read the digit tree, and pass each digit to the member function
 670   // ProcessDigit.
 671   if (!fArrayA) {
 672     AliError("No raw digit array defined");
 673     return kFALSE;
 674   }
 675
 676   Int_t nDigit = fArrayA->GetEntries();
 677   if (nDigit <= 0) return kTRUE;
 678   for (Int_t j = 0; j < nDigit; j++) {
 679     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 680     if (!digit) continue;
 681     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 682       digit->Print();
 683     if (!ProcessRawDigit(digit)) return kFALSE;
 684   }
 685   return kTRUE;
 686 }
 687
 688 //____________________________________________________________________
 689 Bool_t
 690 AliFMDInput::ProcessRecPoints()
 691 {
 692   // Read the reconstrcted points tree, and pass each reconstruction
 693   // object (AliFMDRecPoint) to either ProcessRecPoint.
 694   if (!fTreeR) {
 695     AliError("No recpoint tree defined");
 696     return kFALSE;
 697   }
 698   if (!fArrayR) {
 699     AliError("No recpoints array defined");
 700     return kFALSE;
 701   }
 702
 703   Int_t nEv = fTreeR->GetEntries();
 704   for (Int_t i = 0; i < nEv; i++) {
 705     Int_t recRead  = fTreeR->GetEntry(i);
 706     if (recRead <= 0) continue;
 707     Int_t nRecPoint = fArrayR->GetEntries();
 708     for (Int_t j = 0; j < nRecPoint; j++) {
 709       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 710       if (!recPoint) continue;
 711       if (!ProcessRecPoint(recPoint)) return kFALSE;
 712     }
 713   }
 714   return kTRUE;
 715 }
 716
 717 //____________________________________________________________________
 718 Bool_t
 719 AliFMDInput::ProcessESDs()
 720 {
 721   // Process event summary data
 722   if (!fESD) return kFALSE;
 723   for (UShort_t det = 1; det <= 3; det++) {
 724     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 725     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 726       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 727       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 728       for (UShort_t sec = 0; sec < nsec; sec++) {
 729         for (UShort_t str = 0; str < nstr; str++) {
 730           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 731           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 732           if (!fESD->IsAngleCorrected())
 733             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 734           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 735         }
 736       }
 737     }
 738   }
 739   return kTRUE;
 740 }
 741
 742 //____________________________________________________________________
 743 Bool_t
 744 AliFMDInput::End()
 745 {
 746   // Called at the end of each event.  Per default, it unloads the
 747   // data trees and resets the pointers to the output arrays.   Users
 748   // can overload this, but should call this member function in the
 749   // overloaded member function of the derived class.
 750
 751   // Check if we have been initialized
 752   if (!fIsInit) {
 753     AliError("Not initialized");
 754     return fIsInit;
 755   }
 756   // Possibly unload global kinematics information
 757   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
 758     fLoader->UnloadKinematics();
 759     // fTreeK = 0;
 760     fStack = 0;
 761   }
 762   // Possibly unload FMD Hit information
 763   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
 764     fFMDLoader->UnloadHits();
 765     fTreeH = 0;
 766   }
 767   // Possibly unload FMD Digit information
 768   if (TESTBIT(fTreeMask, kDigits)) {
 769     fFMDLoader->UnloadDigits();
 770     fTreeD = 0;
 771   }
 772   // Possibly unload FMD Sdigit information
 773   if (TESTBIT(fTreeMask, kSDigits)) {
 774     fFMDLoader->UnloadSDigits();
 775     fTreeS = 0;
 776   }
 777   // Possibly unload FMD RecPoints information
 778   if (TESTBIT(fTreeMask, kRecPoints)) {
 779     fFMDLoader->UnloadRecPoints();
 780     fTreeR = 0;
 781   }
 782   // AliInfo("Now out event");
 783   return kTRUE;
 784 }
 785
 786 //____________________________________________________________________
 787 Bool_t
 788 AliFMDInput::Run()
 789 {
 790   // Run over all events and files references in galice.root
 791
 792   Bool_t retval;
 793   if (!(retval = Init())) return retval;
 794
 795   Int_t nEvents = NEvents();
 796   for (Int_t event = 0; nEvents < 0 || event < nEvents; event++) {
 797     if (!(retval = Begin(event))) break;
 798     if (!(retval = Event())) break;
 799     if (!(retval = End())) break;
 800   }
 801   if (!retval) return retval;
 802   retval = Finish();
 803   return retval;
 804 }
 805
 806 //__________________________________________________________________
 807 TArrayF
 808 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
 809 {
 810   // Service function to define a logarithmic axis.
 811   // Parameters:
 812   //   n    Number of bins
 813   //   min  Minimum of axis
 814   //   max  Maximum of axis
 815   TArrayF bins(n+1);
 816   bins[0]      = min;
 817   if (n <= 20) {
 818     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
 819     return bins;
 820   }
 821   Float_t dp   = n / TMath::Log10(max / min);
 822   Float_t pmin = TMath::Log10(min);
 823   for (Int_t i = 1; i < n+1; i++) {
 824     Float_t p = pmin + i / dp;
 825     bins[i]   = TMath::Power(10, p);
 826   }
 827   return bins;
 828 }
 829
 830
 831
 832 //____________________________________________________________________
 833 //
 834 // EOF
 835 //