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