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     fNEvents(-1)
 110 {
 111
 112   // Constructor of an FMD input object.  Specify what data to read in
 113   // using the AddLoad member function.   Sub-classes should at a
 114   // minimum overload the member function Event.   A full job can be
 115   // executed using the member function Run.
 116 }
 117
 118
 119
 120 //____________________________________________________________________
 121 AliFMDInput::AliFMDInput(const char* gAliceFile)
 122   : TNamed("AliFMDInput", "Input handler for various FMD data"),
 123     fGAliceFile(gAliceFile),
 124     fLoader(0),
 125     fRun(0),
 126     fStack(0),
 127     fFMDLoader(0),
 128     fReader(0),
 129     fFMDReader(0),
 130     fFMD(0),
 131     fESD(0),
 132     fESDEvent(0),
 133     fTreeE(0),
 134     fTreeH(0),
 135     fTreeTR(0),
 136     fTreeD(0),
 137     fTreeS(0),
 138     fTreeR(0),
 139     fTreeA(0),
 140     fChainE(0),
 141     fArrayE(0),
 142     fArrayH(0),
 143     fArrayTR(0),
 144     fArrayD(0),
 145     fArrayS(0),
 146     fArrayR(0),
 147     fArrayA(0),
 148     fHeader(0),
 149     fGeoManager(0),
 150     fTreeMask(0),
 151     fRawFile(""),
 152     fIsInit(kFALSE),
 153     fEventCount(0),
 154     fNEvents(-1)
 155 {
 156
 157   // Constructor of an FMD input object.  Specify what data to read in
 158   // using the AddLoad member function.   Sub-classes should at a
 159   // minimum overload the member function Event.   A full job can be
 160   // executed using the member function Run.
 161 }
 162
 163 //____________________________________________________________________
 164 Int_t
 165 AliFMDInput::NEvents() const
 166 {
 167   // Get number of events
 168   if (TESTBIT(fTreeMask, kRaw) ||
 169       TESTBIT(fTreeMask, kRawCalib)) return fReader->GetNumberOfEvents();
 170   if (fChainE) return fChainE->GetEntriesFast();
 171   if (fTreeE) return fTreeE->GetEntries();
 172   return -1;
 173 }
 174
 175 //____________________________________________________________________
 176 Bool_t
 177 AliFMDInput::Init()
 178 {
 179   // Initialize the object.  Get the needed loaders, and such.
 180
 181   // Check if we have been initialized
 182   if (fIsInit) {
 183     AliWarning("Already initialized");
 184     return fIsInit;
 185   }
 186   Info("Init","Initialising w/mask 0x%04x\n"
 187        "\tHits:          %d\n"
 188        "\tKinematics:    %d\n"
 189        "\tDigits:        %d\n"
 190        "\tSDigits:       %d\n"
 191        "\tHeader:        %d\n"
 192        "\tRecPoints:     %d\n"
 193        "\tESD:           %d\n"
 194        "\tRaw:           %d\n"
 195        "\tRawCalib:      %d\n"
 196        "\tGeometry:      %d\n"
 197        "\tTracksRefs:    %d",
 198        fTreeMask,
 199        TESTBIT(fTreeMask, kHits),
 200        TESTBIT(fTreeMask, kKinematics),
 201        TESTBIT(fTreeMask, kDigits),
 202        TESTBIT(fTreeMask, kSDigits),
 203        TESTBIT(fTreeMask, kHeader),
 204        TESTBIT(fTreeMask, kRecPoints),
 205        TESTBIT(fTreeMask, kESD),
 206        TESTBIT(fTreeMask, kRaw),
 207        TESTBIT(fTreeMask, kRawCalib),
 208        TESTBIT(fTreeMask, kGeometry),
 209        TESTBIT(fTreeMask, kTrackRefs));
 210   // Get the run
 211   if (TESTBIT(fTreeMask, kDigits)     ||
 212       TESTBIT(fTreeMask, kSDigits)    ||
 213       TESTBIT(fTreeMask, kKinematics) ||
 214       TESTBIT(fTreeMask, kTrackRefs)  ||
 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     TString fname;
 297     if (fRun) {
 298       fname = gSystem->DirName(fGAliceFile);
 299       fname.Append("/geometry.root");
 300     }
 301     if (!gSystem->AccessPathName(fname.Data()))
 302       fname = "";
 303     AliCDBManager* cdb   = AliCDBManager::Instance();
 304     if (!cdb->IsDefaultStorageSet()) {
 305       cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 306       cdb->SetRun(0);
 307     }
 308
 309     AliGeomManager::LoadGeometry(fname.IsNull() ? 0 : fname.Data());
 310
 311     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
 312     if (align) {
 313       AliInfo("Got alignment data from CDB");
 314       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
 315       if (!array) {
 316         AliWarning("Invalid align data from CDB");
 317       }
 318       else {
 319         Int_t nAlign = array->GetEntries();
 320         for (Int_t i = 0; i < nAlign; i++) {
 321           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
 322           if (!a->ApplyToGeometry()) {
 323             AliWarning(Form("Failed to apply alignment to %s",
 324                             a->GetSymName()));
 325           }
 326         }
 327       }
 328     }
 329     AliFMDGeometry* geom = AliFMDGeometry::Instance();
 330     geom->Init();
 331     geom->InitTransformations();
 332   }
 333
 334   fEventCount = 0;
 335   fIsInit = kTRUE;
 336   return fIsInit;
 337 }
 338
 339 //____________________________________________________________________
 340 Bool_t
 341 AliFMDInput::Begin(Int_t event)
 342 {
 343   // Called at the begining of each event.  Per default, it gets the
 344   // data trees and gets pointers to the output arrays.   Users can
 345   // overload this, but should call this member function in the
 346   // overloaded member function of the derived class.
 347
 348   // Check if we have been initialized
 349   if (!fIsInit) {
 350     AliError("Not initialized");
 351     return fIsInit;
 352   }
 353
 354   // Get the event
 355   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
 356
 357   // Possibly load global kinematics information
 358   if (TESTBIT(fTreeMask, kKinematics)) {
 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)) {
 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)) {
 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   }
 430
 431   // Possibly load FMD Digit information
 432   if (TESTBIT(fTreeMask, kRaw) || TESTBIT(fTreeMask, kRawCalib)) {
 433     Bool_t mon = fRawFile.Contains("mem://");
 434     // AliInfo("Getting FMD raw data digits");
 435     if (mon) std::cout << "Waiting for event ..." << std::flush;
 436     do {
 437       if (!fReader->NextEvent()) {
 438         if (mon) {
 439           gSystem->Sleep(3);
 440           continue;
 441         }
 442         return kFALSE;
 443       }
 444       UInt_t eventType = fReader->GetType();
 445       if(eventType == AliRawEventHeaderBase::kPhysicsEvent ||
 446          eventType == AliRawEventHeaderBase::kCalibrationEvent)
 447         break;
 448     } while (true);
 449     if (mon) std::cout << "got it" << std::endl;
 450     // AliFMDRawReader r(fReader, 0);
 451     fArrayA->Clear();
 452     fFMDReader->ReadAdcs(fArrayA);
 453     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
 454   }
 455   fEventCount++;
 456   return kTRUE;
 457 }
 458
 459
 460 //____________________________________________________________________
 461 Bool_t
 462 AliFMDInput::Event()
 463 {
 464   // Process one event.  The default implementation one or more of
 465   //
 466   //   -  ProcessHits       if the hits are loaded.
 467   //   -  ProcessDigits     if the digits are loaded.
 468   //   -  ProcessSDigits    if the sumbable digits are loaded.
 469   //   -  ProcessRecPoints  if the reconstructed points are loaded.
 470   //   -  ProcessESD        if the event summary data is loaded
 471   //
 472   if (TESTBIT(fTreeMask, kHits))     if (!ProcessHits())      return kFALSE;
 473   if (TESTBIT(fTreeMask, kTrackRefs))if (!ProcessTrackRefs()) return kFALSE;
 474   if (TESTBIT(fTreeMask, kKinematics) &&
 475       TESTBIT(fTreeMask, kHits))     if (!ProcessTracks())    return kFALSE;
 476   if (TESTBIT(fTreeMask, kKinematics))if (!ProcessStack())     return kFALSE;
 477   if (TESTBIT(fTreeMask, kSDigits))  if (!ProcessSDigits())   return kFALSE;
 478   if (TESTBIT(fTreeMask, kDigits))   if (!ProcessDigits())    return kFALSE;
 479   if (TESTBIT(fTreeMask, kRaw))      if (!ProcessRawDigits()) return kFALSE;
 480   if (TESTBIT(fTreeMask, kRawCalib)) if (!ProcessRawCalibDigits())return kFALSE;
 481   if (TESTBIT(fTreeMask, kRecPoints))if (!ProcessRecPoints()) return kFALSE;
 482   if (TESTBIT(fTreeMask, kESD))      if (!ProcessESDs())      return kFALSE;
 483   if (TESTBIT(fTreeMask, kUser))     if (!ProcessUsers())     return kFALSE;
 484
 485   return kTRUE;
 486 }
 487
 488 //____________________________________________________________________
 489 Bool_t
 490 AliFMDInput::ProcessHits()
 491 {
 492   // Read the hit tree, and pass each hit to the member function
 493   // ProcessHit.
 494   if (!fTreeH) {
 495     AliError("No hit tree defined");
 496     return kFALSE;
 497   }
 498   if (!fArrayH) {
 499     AliError("No hit array defined");
 500     return kFALSE;
 501   }
 502
 503   Int_t nTracks = fTreeH->GetEntries();
 504   for (Int_t i = 0; i < nTracks; i++) {
 505     Int_t hitRead  = fTreeH->GetEntry(i);
 506     if (hitRead <= 0) continue;
 507
 508     Int_t nHit = fArrayH->GetEntries();
 509     if (nHit <= 0) continue;
 510
 511     for (Int_t j = 0; j < nHit; j++) {
 512       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 513       if (!hit) continue;
 514
 515       TParticle* track = 0;
 516       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 517         Int_t trackno = hit->Track();
 518         track = fStack->Particle(trackno);
 519       }
 520       if (!ProcessHit(hit, track)) return kFALSE;
 521     }
 522   }
 523   return kTRUE;
 524 }
 525 //____________________________________________________________________
 526 Bool_t
 527 AliFMDInput::ProcessTrackRefs()
 528 {
 529   // Read the reconstrcted points tree, and pass each reconstruction
 530   // object (AliFMDRecPoint) to either ProcessRecPoint.
 531   if (!fTreeTR) {
 532     AliError("No track reference tree defined");
 533     return kFALSE;
 534   }
 535   if (!fArrayTR) {
 536     AliError("No track reference array defined");
 537     return kFALSE;
 538   }
 539
 540   Int_t nEv = fTreeTR->GetEntries();
 541   for (Int_t i = 0; i < nEv; i++) {
 542     Int_t trRead  = fTreeTR->GetEntry(i);
 543     if (trRead <= 0) continue;
 544     Int_t nTrackRefs = fArrayTR->GetEntries();
 545     for (Int_t j = 0; j < nTrackRefs; j++) {
 546       AliTrackReference* trackRef =
 547         static_cast<AliTrackReference*>(fArrayTR->At(j));
 548       if (!trackRef) continue;
 549       // if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
 550       TParticle* track = 0;
 551       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 552         Int_t trackno = trackRef->GetTrack();
 553         track = fStack->Particle(trackno);
 554       }
 555       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
 556     }
 557   }
 558   return kTRUE;
 559 }
 560 //____________________________________________________________________
 561 Bool_t
 562 AliFMDInput::ProcessTracks()
 563 {
 564   // Read the hit tree, and pass each hit to the member function
 565   // ProcessTrack.
 566   if (!fStack) {
 567     AliError("No track tree defined");
 568     return kFALSE;
 569   }
 570   if (!fTreeH) {
 571     AliError("No hit tree defined");
 572     return kFALSE;
 573   }
 574   if (!fArrayH) {
 575     AliError("No hit array defined");
 576     return kFALSE;
 577   }
 578
 579   // Int_t nTracks = fStack->GetNtrack();
 580   Int_t nTracks = fTreeH->GetEntries();
 581   for (Int_t i = 0; i < nTracks; i++) {
 582     Int_t      trackno = nTracks - i - 1;
 583     TParticle* track   = fStack->Particle(trackno);
 584     if (!track) continue;
 585
 586     // Get the hits for this track.
 587     Int_t hitRead  = fTreeH->GetEntry(i);
 588     Int_t nHit     = fArrayH->GetEntries();
 589     if (nHit == 0 || hitRead <= 0) {
 590       // Let user code see the track, even if there's no hits.
 591       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
 592       continue;
 593     }
 594
 595     // Loop over the hits corresponding to this track.
 596     for (Int_t j = 0; j < nHit; j++) {
 597       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 598       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
 599     }
 600   }
 601   return kTRUE;
 602 }
 603 //____________________________________________________________________
 604 Bool_t
 605 AliFMDInput::ProcessStack()
 606 {
 607   // Read the hit tree, and pass each hit to the member function
 608   // ProcessTrack.
 609   if (!fStack) {
 610     AliError("No track tree defined");
 611     return kFALSE;
 612   }
 613   Int_t nTracks = fStack->GetNtrack();
 614   for (Int_t i = 0; i < nTracks; i++) {
 615     Int_t      trackno = nTracks - i - 1;
 616     TParticle* track   = fStack->Particle(trackno);
 617     if (!track) continue;
 618
 619     if (!ProcessParticle(trackno, track)) return kFALSE;
 620   }
 621   return kTRUE;
 622 }
 623 //____________________________________________________________________
 624 Bool_t
 625 AliFMDInput::ProcessDigits()
 626 {
 627   // Read the digit tree, and pass each digit to the member function
 628   // ProcessDigit.
 629   if (!fTreeD) {
 630     AliError("No digit tree defined");
 631     return kFALSE;
 632   }
 633   if (!fArrayD) {
 634     AliError("No digit array defined");
 635     return kFALSE;
 636   }
 637
 638   Int_t nEv = fTreeD->GetEntries();
 639   for (Int_t i = 0; i < nEv; i++) {
 640     Int_t digitRead  = fTreeD->GetEntry(i);
 641     if (digitRead <= 0) continue;
 642     Int_t nDigit = fArrayD->GetEntries();
 643     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
 644     if (nDigit <= 0) continue;
 645     for (Int_t j = 0; j < nDigit; j++) {
 646       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 647       if (!digit) continue;
 648       if (!ProcessDigit(digit)) return kFALSE;
 649     }
 650   }
 651   return kTRUE;
 652 }
 653
 654 //____________________________________________________________________
 655 Bool_t
 656 AliFMDInput::ProcessSDigits()
 657 {
 658   // Read the summable digit tree, and pass each sumable digit to the
 659   // member function ProcessSdigit.
 660   if (!fTreeS) {
 661     AliWarning("No sdigit tree defined");
 662     return kTRUE; // Empty SDigits is fine
 663   }
 664   if (!fArrayS) {
 665     AliWarning("No sdigit array defined");
 666     return kTRUE; // Empty SDigits is fine
 667   }
 668
 669   Int_t nEv = fTreeS->GetEntries();
 670   for (Int_t i = 0; i < nEv; i++) {
 671     Int_t sdigitRead  = fTreeS->GetEntry(i);
 672     if (sdigitRead <= 0) {
 673       AliInfo(Form("Read nothing from tree"));
 674       continue;
 675     }
 676     Int_t nSdigit = fArrayS->GetEntriesFast();
 677     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
 678     AliInfo(Form("Got %5d digits for this event", nSdigit));
 679     if (nSdigit <= 0) continue;
 680     for (Int_t j = 0; j < nSdigit; j++) {
 681       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 682       if (!sdigit) continue;
 683       if (!ProcessSDigit(sdigit)) return kFALSE;
 684     }
 685   }
 686   return kTRUE;
 687 }
 688
 689 //____________________________________________________________________
 690 Bool_t
 691 AliFMDInput::ProcessRawDigits()
 692 {
 693   // Read the digit tree, and pass each digit to the member function
 694   // ProcessDigit.
 695   if (!fArrayA) {
 696     AliError("No raw digit array defined");
 697     return kFALSE;
 698   }
 699
 700   Int_t nDigit = fArrayA->GetEntries();
 701   if (nDigit <= 0) return kTRUE;
 702   for (Int_t j = 0; j < nDigit; j++) {
 703     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 704     if (!digit) continue;
 705     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 706       digit->Print();
 707     if (!ProcessRawDigit(digit)) return kFALSE;
 708   }
 709   return kTRUE;
 710 }
 711
 712 //____________________________________________________________________
 713 Bool_t
 714 AliFMDInput::ProcessRawCalibDigits()
 715 {
 716   // Read the digit tree, and pass each digit to the member function
 717   // ProcessDigit.
 718   if (!fArrayA) {
 719     AliError("No raw digit array defined");
 720     return kFALSE;
 721   }
 722
 723   Int_t nDigit = fArrayA->GetEntries();
 724   if (nDigit <= 0) return kTRUE;
 725   for (Int_t j = 0; j < nDigit; j++) {
 726     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 727     if (!digit) continue;
 728     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
 729       digit->Print();
 730     if (!ProcessRawCalibDigit(digit)) return kFALSE;
 731   }
 732   return kTRUE;
 733 }
 734
 735 //____________________________________________________________________
 736 Bool_t
 737 AliFMDInput::ProcessRecPoints()
 738 {
 739   // Read the reconstrcted points tree, and pass each reconstruction
 740   // object (AliFMDRecPoint) to either ProcessRecPoint.
 741   if (!fTreeR) {
 742     AliError("No recpoint tree defined");
 743     return kFALSE;
 744   }
 745   if (!fArrayR) {
 746     AliError("No recpoints array defined");
 747     return kFALSE;
 748   }
 749
 750   Int_t nEv = fTreeR->GetEntries();
 751   for (Int_t i = 0; i < nEv; i++) {
 752     Int_t recRead  = fTreeR->GetEntry(i);
 753     if (recRead <= 0) continue;
 754     Int_t nRecPoint = fArrayR->GetEntries();
 755     for (Int_t j = 0; j < nRecPoint; j++) {
 756       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 757       if (!recPoint) continue;
 758       if (!ProcessRecPoint(recPoint)) return kFALSE;
 759     }
 760   }
 761   return kTRUE;
 762 }
 763
 764 //____________________________________________________________________
 765 Bool_t
 766 AliFMDInput::ProcessESDs()
 767 {
 768   // Process event summary data
 769   if (!fESD) return kFALSE;
 770   for (UShort_t det = 1; det <= 3; det++) {
 771     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 772     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 773       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 774       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 775       for (UShort_t sec = 0; sec < nsec; sec++) {
 776         for (UShort_t str = 0; str < nstr; str++) {
 777           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 778           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 779           if (!fESD->IsAngleCorrected())
 780             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 781           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 782         }
 783       }
 784     }
 785   }
 786   return kTRUE;
 787 }
 788
 789 //____________________________________________________________________
 790 Bool_t
 791 AliFMDInput::ProcessUsers()
 792 {
 793   // Process event summary data
 794   for (UShort_t det = 1; det <= 3; det++) {
 795     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 796     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 797       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 798       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 799       for (UShort_t sec = 0; sec < nsec; sec++) {
 800         for (UShort_t str = 0; str < nstr; str++) {
 801           Float_t v  = GetSignal(det,*rng,sec,str);
 802           if (!ProcessUser(det, *rng, sec, str, v)) continue;
 803         }
 804       }
 805     }
 806   }
 807   return kTRUE;
 808 }
 809
 810 //____________________________________________________________________
 811 Bool_t
 812 AliFMDInput::End()
 813 {
 814   // Called at the end of each event.  Per default, it unloads the
 815   // data trees and resets the pointers to the output arrays.   Users
 816   // can overload this, but should call this member function in the
 817   // overloaded member function of the derived class.
 818
 819   // Check if we have been initialized
 820   if (!fIsInit) {
 821     AliError("Not initialized");
 822     return fIsInit;
 823   }
 824   // Possibly unload global kinematics information
 825   if (TESTBIT(fTreeMask, kKinematics)) {
 826     fLoader->UnloadKinematics();
 827     // fTreeK = 0;
 828     fStack = 0;
 829   }
 830   // Possibly unload FMD Hit information
 831   if (TESTBIT(fTreeMask, kHits)) {
 832     fFMDLoader->UnloadHits();
 833     fTreeH = 0;
 834   }
 835   // Possibly unload FMD Digit information
 836   if (TESTBIT(fTreeMask, kDigits)) {
 837     fFMDLoader->UnloadDigits();
 838     fTreeD = 0;
 839   }
 840   // Possibly unload FMD Sdigit information
 841   if (TESTBIT(fTreeMask, kSDigits)) {
 842     fFMDLoader->UnloadSDigits();
 843     fTreeS = 0;
 844   }
 845   // Possibly unload FMD RecPoints information
 846   if (TESTBIT(fTreeMask, kRecPoints)) {
 847     fFMDLoader->UnloadRecPoints();
 848     fTreeR = 0;
 849   }
 850   // AliInfo("Now out event");
 851   return kTRUE;
 852 }
 853
 854 //____________________________________________________________________
 855 Bool_t
 856 AliFMDInput::Run(UInt_t maxEvents)
 857 {
 858   // Run over all events and files references in galice.root
 859
 860   Bool_t retval;
 861   if (!(retval = Init())) return retval;
 862
 863   fNEvents = NEvents();
 864   if (fNEvents < 0)       fNEvents = maxEvents;
 865   else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
 866
 867   Int_t event = 0;
 868   for (; fNEvents < 0 || event < fNEvents; event++) {
 869     printf("\rEvent %8d/%8d ...", event, fNEvents);
 870     if (!(retval = Begin(event))) break;
 871     if (!(retval = Event())) break;
 872     if (!(retval = End())) break;
 873   }
 874   printf("Looped over %8d events\n", event+1);
 875   if (!retval) return retval;
 876   retval = Finish();
 877   return retval;
 878 }
 879
 880 //__________________________________________________________________
 881 TArrayF
 882 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
 883 {
 884   // Service function to define a logarithmic axis.
 885   // Parameters:
 886   //   n    Number of bins
 887   //   min  Minimum of axis
 888   //   max  Maximum of axis
 889   TArrayF bins(n+1);
 890   bins[0]      = min;
 891   if (n <= 20) {
 892     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
 893     return bins;
 894   }
 895   Float_t dp   = n / TMath::Log10(max / min);
 896   Float_t pmin = TMath::Log10(min);
 897   for (Int_t i = 1; i < n+1; i++) {
 898     Float_t p = pmin + i / dp;
 899     bins[i]   = TMath::Power(10, p);
 900   }
 901   return bins;
 902 }
 903
 904
 905
 906 //____________________________________________________________________
 907 //
 908 // EOF
 909 //