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 "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 <AliAlignObjParams.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           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(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) || TESTBIT(fTreeMask, kTracks)) {
 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) || TESTBIT(fTreeMask, kTracks)) {
 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, kTracks))
 356     if (!ProcessTracks()) return kFALSE;
 357   if (TESTBIT(fTreeMask, kDigits))
 358     if (!ProcessDigits()) return kFALSE;
 359   if (TESTBIT(fTreeMask, kSDigits))
 360     if (!ProcessSDigits()) return kFALSE;
 361   if (TESTBIT(fTreeMask, kRaw))
 362     if (!ProcessRawDigits()) return kFALSE;
 363   if (TESTBIT(fTreeMask, kRecPoints))
 364     if (!ProcessRecPoints()) return kFALSE;
 365   if (TESTBIT(fTreeMask, kESD))
 366     if (!ProcessESDs()) return kFALSE;
 367
 368   return kTRUE;
 369 }
 370
 371 //____________________________________________________________________
 372 Bool_t
 373 AliFMDInput::ProcessHits()
 374 {
 375   // Read the hit tree, and pass each hit to the member function
 376   // ProcessHit.
 377   if (!fTreeH) {
 378     AliError("No hit tree defined");
 379     return kFALSE;
 380   }
 381   Int_t nTracks = fTreeH->GetEntries();
 382   for (Int_t i = 0; i < nTracks; i++) {
 383     Int_t hitRead  = fTreeH->GetEntry(i);
 384     if (hitRead <= 0) continue;
 385     if (!fArrayH) {
 386       AliError("No hit array defined");
 387       return kFALSE;
 388     }
 389     Int_t nHit = fArrayH->GetEntries();
 390     if (nHit <= 0) continue;
 391     for (Int_t j = 0; j < nHit; j++) {
 392       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 393       if (!hit) continue;
 394       TParticle* track = 0;
 395       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
 396         Int_t trackno = hit->Track();
 397         track = fStack->Particle(trackno);
 398       }
 399       if (!ProcessHit(hit, track)) return kFALSE;
 400     }
 401   }
 402   return kTRUE;
 403 }
 404
 405 //____________________________________________________________________
 406 Bool_t
 407 AliFMDInput::ProcessTracks()
 408 {
 409   // Read the hit tree, and pass each hit to the member function
 410   // ProcessHit.
 411   if (!fStack) {
 412     AliError("No track tree defined");
 413     return kFALSE;
 414   }
 415   if (!fTreeH) {
 416     AliError("No hit tree defined");
 417     return kFALSE;
 418   }
 419   Int_t nTracks = fTreeH->GetEntries();
 420   for (Int_t i = 0; i < nTracks; i++) {
 421     TParticle* track = fStack->Particle(i);
 422     if (!track) continue;
 423     Int_t hitRead  = fTreeH->GetEntry(i);
 424     if (hitRead <= 0) continue;
 425     if (!fArrayH) {
 426       AliError("No hit array defined");
 427       return kFALSE;
 428     }
 429     Int_t nHit = fArrayH->GetEntries();
 430     if (nHit <= 0) continue;
 431
 432     for (Int_t j = 0; j < nHit; j++) {
 433       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
 434       if (!hit) continue;
 435       if (!ProcessTrack(i, track, hit)) return kFALSE;
 436     }
 437     // if (!ProcessTrack(i, track, fArrayH)) return kFALSE;
 438   }
 439   return kTRUE;
 440 }
 441
 442 //____________________________________________________________________
 443 Bool_t
 444 AliFMDInput::ProcessDigits()
 445 {
 446   // Read the digit tree, and pass each digit to the member function
 447   // ProcessDigit.
 448   Int_t nEv = fTreeD->GetEntries();
 449   for (Int_t i = 0; i < nEv; i++) {
 450     Int_t digitRead  = fTreeD->GetEntry(i);
 451     if (digitRead <= 0) continue;
 452     Int_t nDigit = fArrayD->GetEntries();
 453     if (nDigit <= 0) continue;
 454     for (Int_t j = 0; j < nDigit; j++) {
 455       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
 456       if (!digit) continue;
 457       if (!ProcessDigit(digit)) return kFALSE;
 458     }
 459   }
 460   return kTRUE;
 461 }
 462
 463 //____________________________________________________________________
 464 Bool_t
 465 AliFMDInput::ProcessSDigits()
 466 {
 467   // Read the summable digit tree, and pass each sumable digit to the
 468   // member function ProcessSdigit.
 469   Int_t nEv = fTreeD->GetEntries();
 470   for (Int_t i = 0; i < nEv; i++) {
 471     Int_t sdigitRead  = fTreeS->GetEntry(i);
 472     if (sdigitRead <= 0) continue;
 473     Int_t nSdigit = fArrayS->GetEntries();
 474     if (nSdigit <= 0) continue;
 475     for (Int_t j = 0; j < nSdigit; j++) {
 476       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
 477       if (!sdigit) continue;
 478       if (!ProcessSDigit(sdigit)) return kFALSE;
 479     }
 480   }
 481   return kTRUE;
 482 }
 483
 484 //____________________________________________________________________
 485 Bool_t
 486 AliFMDInput::ProcessRawDigits()
 487 {
 488   // Read the digit tree, and pass each digit to the member function
 489   // ProcessDigit.
 490   Int_t nDigit = fArrayA->GetEntries();
 491   if (nDigit <= 0) return kTRUE;
 492   for (Int_t j = 0; j < nDigit; j++) {
 493     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
 494     if (!digit) continue;
 495     if (!ProcessRawDigit(digit)) return kFALSE;
 496   }
 497   return kTRUE;
 498 }
 499
 500 //____________________________________________________________________
 501 Bool_t
 502 AliFMDInput::ProcessRecPoints()
 503 {
 504   // Read the reconstrcted points tree, and pass each reconstruction
 505   // object (AliFMDRecPoint) to either ProcessRecPoint.
 506   Int_t nEv = fTreeR->GetEntries();
 507   for (Int_t i = 0; i < nEv; i++) {
 508     Int_t recRead  = fTreeR->GetEntry(i);
 509     if (recRead <= 0) continue;
 510     Int_t nRecPoint = fArrayR->GetEntries();
 511     for (Int_t j = 0; j < nRecPoint; j++) {
 512       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
 513       if (!recPoint) continue;
 514       if (!ProcessRecPoint(recPoint)) return kFALSE;
 515     }
 516   }
 517   return kTRUE;
 518 }
 519
 520 //____________________________________________________________________
 521 Bool_t
 522 AliFMDInput::ProcessESDs()
 523 {
 524   // Process event summary data
 525   if (!fESD) return kFALSE;
 526   for (UShort_t det = 1; det <= 3; det++) {
 527     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
 528     for (Char_t* rng = rings; *rng != '\0'; rng++) {
 529       UShort_t nsec = (*rng == 'I' ?  20 :  40);
 530       UShort_t nstr = (*rng == 'I' ? 512 : 256);
 531       for (UShort_t sec = 0; sec < nsec; sec++) {
 532         for (UShort_t str = 0; str < nstr; str++) {
 533           Float_t eta  = fESD->Eta(det,*rng,sec,str);
 534           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
 535           if (!fESD->IsAngleCorrected())
 536             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
 537           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
 538         }
 539       }
 540     }
 541   }
 542   return kTRUE;
 543 }
 544
 545 //____________________________________________________________________
 546 Bool_t
 547 AliFMDInput::End()
 548 {
 549   // Called at the end of each event.  Per default, it unloads the
 550   // data trees and resets the pointers to the output arrays.   Users
 551   // can overload this, but should call this member function in the
 552   // overloaded member function of the derived class.
 553
 554   // Check if we have been initialized
 555   if (!fIsInit) {
 556     AliError("Not initialized");
 557     return fIsInit;
 558   }
 559   // Possibly unload global kinematics information
 560   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
 561     fLoader->UnloadKinematics();
 562     // fTreeK = 0;
 563     fStack = 0;
 564   }
 565   // Possibly unload FMD Hit information
 566   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
 567     fFMDLoader->UnloadHits();
 568     fTreeH = 0;
 569   }
 570   // Possibly unload FMD Digit information
 571   if (TESTBIT(fTreeMask, kDigits)) {
 572     fFMDLoader->UnloadDigits();
 573     fTreeD = 0;
 574   }
 575   // Possibly unload FMD Sdigit information
 576   if (TESTBIT(fTreeMask, kSDigits)) {
 577     fFMDLoader->UnloadSDigits();
 578     fTreeS = 0;
 579   }
 580   // Possibly unload FMD RecPoints information
 581   if (TESTBIT(fTreeMask, kRecPoints)) {
 582     fFMDLoader->UnloadRecPoints();
 583     fTreeR = 0;
 584   }
 585   AliInfo("Now out event");
 586   return kTRUE;
 587 }
 588
 589 //____________________________________________________________________
 590 Bool_t
 591 AliFMDInput::Run()
 592 {
 593   // Run over all events and files references in galice.root
 594
 595   Bool_t retval;
 596   if (!(retval = Init())) return retval;
 597
 598   Int_t nEvents = NEvents();
 599   for (Int_t event = 0; event < nEvents; event++) {
 600     if (!(retval = Begin(event))) break;
 601     if (!(retval = Event())) break;
 602     if (!(retval = End())) break;
 603   }
 604   if (!retval) return retval;
 605   retval = Finish();
 606   return retval;
 607 }
 608
 609
 610 //____________________________________________________________________
 611 //
 612 // EOF
 613 //