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