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 "AliLog.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 "AliFMD.h"             // ALIFMD_H
  39 #include "AliFMDHit.h"          // ALIFMDHIT_H
  40 #include "AliFMDDigit.h"        // ALIFMDDigit_H
  41 #include "AliFMDSDigit.h"       // ALIFMDDigit_H
  42 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
  43 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
  44 #include <AliESD.h>
  45 #include <AliESDFMD.h>
  46 #include <AliCDBManager.h>
  47 #include <AliCDBEntry.h>
  48 #include <AliAlignObjAngles.h>
  49 #include <TTree.h>              // ROOT_TTree
  50 #include <TChain.h>             // ROOT_TChain
  51 #include <TParticle.h>          // ROOT_TParticle
  52 #include <TString.h>            // ROOT_TString
  53 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
  54 #include <TMath.h>              // ROOT_TMath
  55 #include <TGeoManager.h>        // ROOT_TGeoManager
  56 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
  57 #include <Riostream.h>          // ROOT_Riostream
  58 #include <TFile.h>              // ROOT_TFile
  59 #include <TStreamerInfo.h>
  60 //____________________________________________________________________
  61 ClassImp(AliFMDInput)
  62 #if 0
  63   ; // This is here to keep Emacs for indenting the next line
  64 #endif
  65
  66
  67 //____________________________________________________________________
  68 AliFMDInput::AliFMDInput()
  69   : fGAliceFile(""),
  70     fLoader(0),
  71     fRun(0),
  72     fStack(0),
  73     fFMDLoader(0),
  74     fReader(0),
  75     fFMD(0),
  76     fMainESD(0),
  77     fESD(0),
  78     fTreeE(0),
  79     fTreeH(0),
  80     fTreeD(0),
  81     fTreeS(0),
  82     fTreeR(0),
  83     fTreeA(0),
  84     fChainE(0),
  85     fArrayE(0),
  86     fArrayH(0),
  87     fArrayD(0),
  88     fArrayS(0),
  89     fArrayR(0),
  90     fArrayA(0),
  91     fGeoManager(0),
  92     fTreeMask(0),
  93     fIsInit(kFALSE)
  94 {
  95   // Constructor of an FMD input object.  Specify what data to read in
  96   // using the AddLoad member function.   Sub-classes should at a
  97   // minimum overload the member function Event.   A full job can be
  98   // executed using the member function Run.
  99 }
 100
 101
 102
 103 //____________________________________________________________________
 104 AliFMDInput::AliFMDInput(const char* gAliceFile)
 105   : fGAliceFile(gAliceFile),
 106     fLoader(0),
 107     fRun(0),
 108     fStack(0),
 109     fFMDLoader(0),
 110     fReader(0),
 111     fFMD(0),
 112     fMainESD(0),
 113     fESD(0),
 114     fTreeE(0),
 115     fTreeH(0),
 116     fTreeD(0),
 117     fTreeS(0),
 118     fTreeR(0),
 119     fTreeA(0),
 120     fChainE(0),
 121     fArrayE(0),
 122     fArrayH(0),
 123     fArrayD(0),
 124     fArrayS(0),
 125     fArrayR(0),
 126     fArrayA(0),
 127     fGeoManager(0),
 128     fTreeMask(0),
 129     fIsInit(kFALSE)
 130 {
 131   // Constructor of an FMD input object.  Specify what data to read in
 132   // using the AddLoad member function.   Sub-classes should at a
 133   // minimum overload the member function Event.   A full job can be
 134   // executed using the member function Run.
 135 }
 136
 137 //____________________________________________________________________
 138 Int_t
 139 AliFMDInput::NEvents() const
 140 {
 141   // Get number of events
 142   if (fTreeE) return fTreeE->GetEntries();
 143   return -1;
 144 }
 145
 146 //____________________________________________________________________
 147 Bool_t
 148 AliFMDInput::Init()
 149 {
 150   // Initialize the object.  Get the needed loaders, and such.
 151
 152   // Check if we have been initialized
 153   if (fIsInit) {
 154     AliWarning("Already initialized");
 155     return fIsInit;
 156   }
 157   if (fGAliceFile.IsNull()) fGAliceFile = "galice.root";
 158   // Get the loader
 159   fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
 160   if (!fLoader) {
 161     AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
 162     return kFALSE;
 163   }
 164
 165   // Get the run
 166   if  (fLoader->LoadgAlice()) return kFALSE;
 167   fRun = fLoader->GetAliRun();
 168
 169   // Get the FMD
 170   fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
 171   if (!fFMD) {
 172     AliError("Failed to get detector FMD from loader");
 173     return kFALSE;
 174   }
 175
 176   // Get the FMD loader
 177   fFMDLoader = fLoader->GetLoader("FMDLoader");
 178   if (!fFMDLoader) {
 179     AliError("Failed to get detector FMD loader from loader");
 180     return kFALSE;
 181   }
 182   if (fLoader->LoadHeader()) {
 183     AliError("Failed to get event header information from loader");
 184     return kFALSE;
 185   }
 186   fTreeE = fLoader->TreeE();
 187
 188   // Optionally, get the ESD files
 189   if (TESTBIT(fTreeMask, kESD)) {
 190     fChainE = new TChain("esdTree");
 191     TSystemDirectory dir(".",".");
 192     TList*           files = dir.GetListOfFiles();
 193     TSystemFile*     file = 0;
 194     if (!files) {
 195       AliError("No files");
 196       return kFALSE;
 197     }
 198     files->Sort();
 199     TIter            next(files);
 200     while ((file = static_cast<TSystemFile*>(next()))) {
 201       TString fname(file->GetName());
 202       if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
 203     }
 204     fChainE->SetBranchAddress("ESD", &fMainESD);
 205   }
 206
 207   if (TESTBIT(fTreeMask, kRaw)) {
 208     AliInfo("Getting FMD raw data digits");
 209     fArrayA = new TClonesArray("AliFMDDigit");
 210     fReader = new AliRawReaderFile(-1);
 211   }
 212
 213   // Optionally, get the geometry
 214   if (TESTBIT(fTreeMask, kGeometry)) {
 215     TString fname(fRun->GetGeometryFileName());
 216     if (fname.IsNull()) {
 217       Warning("Init", "No file name for the geometry from AliRun");
 218       fname = gSystem->DirName(fGAliceFile);
 219       fname.Append("/geometry.root");
 220     }
 221     fGeoManager = TGeoManager::Import(fname.Data());
 222     if (!fGeoManager) {
 223       Fatal("Init", "No geometry manager found");
 224       return kFALSE;
 225     }
 226     AliCDBManager* cdb   = AliCDBManager::Instance();
 227     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
 228     if (align) {
 229       AliInfo("Got alignment data from CDB");
 230       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
 231       if (!array) {
 232         AliWarning("Invalid align data from CDB");
 233       }
 234       else {
 235         Int_t nAlign = array->GetEntries();
 236         for (Int_t i = 0; i < nAlign; i++) {
 237           AliAlignObjAngles* a = static_cast<AliAlignObjAngles*>(array->At(i));
 238           if (!a->ApplyToGeometry()) {
 239             AliWarning(Form("Failed to apply alignment to %s",
 240                             a->GetSymName()));
 241           }
 242         }
 243       }
 244     }
 245   }
 246
 247
 248   fIsInit = kTRUE;
 249   return fIsInit;
 250 }
 251
 252 //____________________________________________________________________
 253 Bool_t
 254 AliFMDInput::Begin(Int_t event)
 255 {
 256   // Called at the begining of each event.  Per default, it gets the
 257   // data trees and gets pointers to the output arrays.   Users can
 258   // overload this, but should call this member function in the
 259   // overloaded member function of the derived class.
 260
 261   // Check if we have been initialized
 262   if (!fIsInit) {
 263     AliError("Not initialized");
 264     return fIsInit;
 265   }
 266   // Get the event
 267   if (fLoader->GetEvent(event)) return kFALSE;
 268   AliInfo(Form("Now in event %d/%d", event, NEvents()));
 269
 270   // Possibly load global kinematics information
 271   if (TESTBIT(fTreeMask, kKinematics)) {
 272     AliInfo("Getting kinematics");
 273     if (fLoader->LoadKinematics()) return kFALSE;
 274     fStack = fLoader->Stack();
 275   }
 276   // Possibly load FMD Hit information
 277   if (TESTBIT(fTreeMask, kHits)) {
 278     AliInfo("Getting FMD hits");
 279     if (fFMDLoader->LoadHits()) return kFALSE;
 280     fTreeH = fFMDLoader->TreeH();
 281     if (!fArrayH) fArrayH = fFMD->Hits();
 282   }
 283   // Possibly load FMD Digit information
 284   if (TESTBIT(fTreeMask, kDigits)) {
 285     AliInfo("Getting FMD digits");
 286     if (fFMDLoader->LoadDigits()) return kFALSE;
 287     fTreeD = fFMDLoader->TreeD();
 288     if (fTreeD) {
 289       if (!fArrayD) fArrayD = fFMD->Digits();
 290     }
 291     else {
 292       fArrayD = 0;
 293       AliWarning(Form("Failed to load FMD Digits"));
 294     }
 295   }
 296   // Possibly load FMD Sdigit information
 297   if (TESTBIT(fTreeMask, kSDigits)) {
 298     AliInfo("Getting FMD summable digits");
 299     if (fFMDLoader->LoadSDigits()) return kFALSE;
 300     fTreeS = fFMDLoader->TreeS();
 301     if (!fArrayS) fArrayS = fFMD->SDigits();
 302   }
 303   // Possibly load FMD RecPoints information
 304   if (TESTBIT(fTreeMask, kRecPoints)) {
 305     AliInfo("Getting FMD reconstructed points");
 306     if (fFMDLoader->LoadRecPoints()) return kFALSE;
 307     fTreeR = fFMDLoader->TreeR();
 308     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
 309     fTreeR->SetBranchAddress("FMD",  &fArrayR);
 310   }  // Possibly load FMD ESD information
 311   if (TESTBIT(fTreeMask, kESD)) {
 312     AliInfo("Getting FMD event summary data");
 313     Int_t read = fChainE->GetEntry(event);
 314     if (read <= 0) return kFALSE;
 315     fESD = fMainESD->GetFMDData();
 316     if (!fESD) return kFALSE;
 317     TFile* f = fChainE->GetFile();
 318     if (f) {
 319       TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
 320       if (o) {
 321         TStreamerInfo* info = static_cast<TStreamerInfo*>(o);
 322         std::cout << "AliFMDMap class version read is "
 323                   <<  info->GetClassVersion() << std::endl;
 324       }
 325     }
 326     // fESD->CheckNeedUShort(fChainE->GetFile());
 327   }
 328   // Possibly load FMD Digit information
 329   if (TESTBIT(fTreeMask, kRaw)) {
 330     AliInfo("Getting FMD raw data digits");
 331     if (!fReader->NextEvent()) return kFALSE;
 332     AliFMDRawReader r(fReader, 0);
 333     fArrayA->Clear();
 334     r.ReadAdcs(fArrayA);
 335   }
 336
 337   return kTRUE;
 338 }
 339
 340
 341 //____________________________________________________________________
 342 Bool_t
 343 AliFMDInput::Event()
 344 {
 345   // Process one event.  The default implementation one or more of
 346   //
 347   //   -  ProcessHits       if the hits are loaded.
 348   //   -  ProcessDigits     if the digits are loaded.
 349   //   -  ProcessSDigits    if the sumbable digits are loaded.
 350   //   -  ProcessRecPoints  if the reconstructed points are loaded.
 351   //   -  ProcessESD        if the event summary data is loaded
 352   //
 353   if (TESTBIT(fTreeMask, kHits))
 354     if (!ProcessHits()) return kFALSE;
 355   if (TESTBIT(fTreeMask, kDigits))
 356     if (!ProcessDigits()) return kFALSE;
 357   if (TESTBIT(fTreeMask, kSDigits))
 358     if (!ProcessSDigits()) return kFALSE;
 359   if (TESTBIT(fTreeMask, kRaw))
 360     if (!ProcessRawDigits()) return kFALSE;
 361   if (TESTBIT(fTreeMask, kRecPoints))
 362     if (!ProcessRecPoints()) return kFALSE;
 363   if (TESTBIT(fTreeMask, kESD))
 364     if (!ProcessESDs()) return kFALSE;
 365
 366   return kTRUE;
 367 }
 368
 369 //____________________________________________________________________
 370 Bool_t
 371 AliFMDInput::ProcessHits()
 372 {
 373   // Read the hit tree, and pass each hit to the member function
 374   // ProcessHit.
 375   if (!fTreeH) {
 376     AliError("No hit tree defined");
 377     return kFALSE;
 378   }
 379   Int_t nTracks = fTreeH->GetEntries();
 380   for (Int_t i = 0; i < nTracks; i++) {
 381     Int_t hitRead  = fTreeH->GetEntry(i);
 382     if (hitRead <= 0) continue;
 383     if (!fArrayH) {
 384       AliError("No hit array defined");
 385       return kFALSE;
 386     }
 387     Int_t nHit = fArrayH->GetEntries();
 388     if (nHit <= 0) continue;
 389     for (Int_t j = 0; j < nHit; j++) {
 390       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 391       if (!hit) continue;
 392       TParticle* track = 0;
 393       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 394         Int_t trackno = hit->Track();
 395         track = fStack->Particle(trackno);
 396       }
 397       if (!ProcessHit(hit, track)) return kFALSE;
 398     }
 399   }
 400   return kTRUE;
 401 }
 402
 403 //____________________________________________________________________
 404 Bool_t
 405 AliFMDInput::ProcessDigits()
 406 {
 407   // Read the digit tree, and pass each digit to the member function
 408   // ProcessDigit.
 409   Int_t nEv = fTreeD->GetEntries();
 410   for (Int_t i = 0; i < nEv; i++) {
 411     Int_t digitRead  = fTreeD->GetEntry(i);
 412     if (digitRead <= 0) continue;
 413     Int_t nDigit = fArrayD->GetEntries();
 414     if (nDigit <= 0) continue;
 415     for (Int_t j = 0; j < nDigit; j++) {
 416       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 417       if (!digit) continue;
 418       if (!ProcessDigit(digit)) return kFALSE;
 419     }
 420   }
 421   return kTRUE;
 422 }
 423
 424 //____________________________________________________________________
 425 Bool_t
 426 AliFMDInput::ProcessSDigits()
 427 {
 428   // Read the summable digit tree, and pass each sumable digit to the
 429   // member function ProcessSdigit.
 430   Int_t nEv = fTreeD->GetEntries();
 431   for (Int_t i = 0; i < nEv; i++) {
 432     Int_t sdigitRead  = fTreeS->GetEntry(i);
 433     if (sdigitRead <= 0) continue;
 434     Int_t nSdigit = fArrayS->GetEntries();
 435     if (nSdigit <= 0) continue;
 436     for (Int_t j = 0; j < nSdigit; j++) {
 437       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 438       if (!sdigit) continue;
 439       if (!ProcessSDigit(sdigit)) return kFALSE;
 440     }
 441   }
 442   return kTRUE;
 443 }
 444
 445 //____________________________________________________________________
 446 Bool_t
 447 AliFMDInput::ProcessRawDigits()
 448 {
 449   // Read the digit tree, and pass each digit to the member function
 450   // ProcessDigit.
 451   Int_t nDigit = fArrayA->GetEntries();
 452   if (nDigit <= 0) return kTRUE;
 453   for (Int_t j = 0; j < nDigit; j++) {
 454     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 455     if (!digit) continue;
 456     if (!ProcessRawDigit(digit)) return kFALSE;
 457   }
 458   return kTRUE;
 459 }
 460
 461 //____________________________________________________________________
 462 Bool_t
 463 AliFMDInput::ProcessRecPoints()
 464 {
 465   // Read the reconstrcted points tree, and pass each reconstruction
 466   // object (AliFMDRecPoint) to either ProcessRecPoint.
 467   Int_t nEv = fTreeR->GetEntries();
 468   for (Int_t i = 0; i < nEv; i++) {
 469     Int_t recRead  = fTreeR->GetEntry(i);
 470     if (recRead <= 0) continue;
 471     Int_t nRecPoint = fArrayR->GetEntries();
 472     for (Int_t j = 0; j < nRecPoint; j++) {
 473       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 474       if (!recPoint) continue;
 475       if (!ProcessRecPoint(recPoint)) return kFALSE;
 476     }
 477   }
 478   return kTRUE;
 479 }
 480
 481 //____________________________________________________________________
 482 Bool_t
 483 AliFMDInput::ProcessESDs()
 484 {
 485   // Process event summary data
 486   if (!fESD) return kFALSE;
 487   for (UShort_t det = 1; det <= 3; det++) {
 488     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 489     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 490       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 491       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 492       for (UShort_t sec = 0; sec < nsec; sec++) {
 493         for (UShort_t str = 0; str < nstr; str++) {
 494           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 495           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 496           if (!fESD->IsAngleCorrected())
 497             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 498           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 499         }
 500       }
 501     }
 502   }
 503   return kTRUE;
 504 }
 505
 506 //____________________________________________________________________
 507 Bool_t
 508 AliFMDInput::End()
 509 {
 510   // Called at the end of each event.  Per default, it unloads the
 511   // data trees and resets the pointers to the output arrays.   Users
 512   // can overload this, but should call this member function in the
 513   // overloaded member function of the derived class.
 514
 515   // Check if we have been initialized
 516   if (!fIsInit) {
 517     AliError("Not initialized");
 518     return fIsInit;
 519   }
 520   // Possibly unload global kinematics information
 521   if (TESTBIT(fTreeMask, kKinematics)) {
 522     fLoader->UnloadKinematics();
 523     // fTreeK = 0;
 524     fStack = 0;
 525   }
 526   // Possibly unload FMD Hit information
 527   if (TESTBIT(fTreeMask, kHits)) {
 528     fFMDLoader->UnloadHits();
 529     fTreeH = 0;
 530   }
 531   // Possibly unload FMD Digit information
 532   if (TESTBIT(fTreeMask, kDigits)) {
 533     fFMDLoader->UnloadDigits();
 534     fTreeD = 0;
 535   }
 536   // Possibly unload FMD Sdigit information
 537   if (TESTBIT(fTreeMask, kSDigits)) {
 538     fFMDLoader->UnloadSDigits();
 539     fTreeS = 0;
 540   }
 541   // Possibly unload FMD RecPoints information
 542   if (TESTBIT(fTreeMask, kRecPoints)) {
 543     fFMDLoader->UnloadRecPoints();
 544     fTreeR = 0;
 545   }
 546   AliInfo("Now out event");
 547   return kTRUE;
 548 }
 549
 550 //____________________________________________________________________
 551 Bool_t
 552 AliFMDInput::Run()
 553 {
 554   // Run over all events and files references in galice.root
 555
 556   Bool_t retval;
 557   if (!(retval = Init())) return retval;
 558
 559   Int_t nEvents = NEvents();
 560   for (Int_t event = 0; event < nEvents; event++) {
 561     if (!(retval = Begin(event))) break;
 562     if (!(retval = Event())) break;
 563     if (!(retval = End())) break;
 564   }
 565   if (!retval) return retval;
 566   retval = Finish();
 567   return retval;
 568 }
 569
 570
 571 //____________________________________________________________________
 572 //
 573 // EOF
 574 //