PHOS/AliPHOSGetter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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
  16 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 //  A singleton. This class should be used in the analysis stage to get
  20 //  reconstructed objects: Digits, RecPoints, TrackSegments and RecParticles,
  21 //  instead of directly reading them from galice.root file. This container
  22 //  ensures, that one reads Digits, made of these particular digits, RecPoints,
  23 //  made of these particular RecPoints, TrackSegments and RecParticles.
  24 //  This becomes non trivial if there are several identical branches, produced with
  25 //  different set of parameters.
  26 //
  27 //  An example of how to use (see also class AliPHOSAnalyser):
  28 //  AliPHOSGetter * gime = AliPHOSGetter::GetInstance("galice.root","test") ;
  29 //  for(Int_t irecp = 0; irecp < gime->NRecParticles() ; irecp++)
  30 //     AliPHOSRecParticle * part = gime->RecParticle(1) ;
  31 //     ................
  32 //  gime->Event(event) ;    // reads new event from galice.root
  33 //
  34 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (RRC KI & SUBATECH)
  35 //*--         Completely redesigned by Dmitri Peressounko March 2001
  36 //
  37 //*-- YS June 2001 : renamed the original AliPHOSIndexToObject and make
  38 //*--         systematic usage of TFolders without changing the interface
  39 //////////////////////////////////////////////////////////////////////////////
  40
  41 // --- ROOT system ---
  42
  43 #include <TFile.h>
  44 #include <TROOT.h>
  45 #include <TSystem.h>
  46 #include <TParticle.h>
  47 #include <TF1.h>
  48 #include <TGraph.h>
  49 #include <TCanvas.h>
  50 #include <TStyle.h>
  51 //#include <TFrame.h>
  52
  53 // --- Standard library ---
  54
  55 // --- AliRoot header files ---
  56 #include "AliESD.h"
  57 #include "AliHeader.h"
  58 #include "AliMC.h"
  59 #include "AliPHOS.h"
  60 #include "AliPHOSBeamTestEvent.h"
  61 #include "AliPHOSGetter.h"
  62 #include "AliPHOSLoader.h"
  63 #include "AliRunLoader.h"
  64 #include "AliStack.h"
  65 #include "AliCaloRawStream.h"
  66 #include "AliRawReaderFile.h"
  67 #include "AliLog.h"
  68 #include "AliCDBLocal.h"
  69 #include "AliCDBStorage.h"
  70 #include "AliCDBManager.h"
  71
  72 ClassImp(AliPHOSGetter)
  73
  74 AliPHOSGetter   * AliPHOSGetter::fgObjGetter  = 0 ;
  75 AliPHOSLoader   * AliPHOSGetter::fgPhosLoader = 0;
  76 AliPHOSCalibData* AliPHOSGetter::fgCalibData  = 0;
  77 Int_t AliPHOSGetter::fgDebug = 0;
  78
  79 //  TFile * AliPHOSGetter::fgFile = 0 ;
  80
  81 AliPHOSGetter::AliPHOSGetter() :
  82   fBTE(0),
  83   fLoadingStatus(),
  84   fNPrimaries(0),
  85   fPrimaries(0),
  86   fESDFile(0),
  87   fESDFileName(),
  88   fESD(0),
  89   fESDTree(0),
  90   fRawDigits(kFALSE),
  91   fcdb(0)
  92 {
  93   // ctor: this is a singleton, the ctor should never be called but cint needs it as public
  94   Fatal("ctor", "AliPHOSGetter is a singleton default ctor not callable") ;
  95 }
  96
  97
  98 //____________________________________________________________________________
  99 AliPHOSGetter::AliPHOSGetter(const char* headerFile, const char* version, Option_t * openingOption) :
 100   fBTE(0),
 101   fLoadingStatus(),
 102   fNPrimaries(0),
 103   fPrimaries(0),
 104   fESDFile(0),
 105   fESDFileName(),
 106   fESD(0),
 107   fESDTree(0),
 108   fRawDigits(kFALSE),
 109   fcdb(0)
 110 {
 111   // ctor only called by Instance()
 112
 113   AliRunLoader* rl = AliRunLoader::GetRunLoader(version) ;
 114   if (!rl) {
 115     rl = AliRunLoader::Open(headerFile, version, openingOption);
 116     if (!rl) {
 117       Fatal("AliPHOSGetter", "Could not find the Run Loader for %s - %s",headerFile, version) ;
 118       return ;
 119     }
 120     if (rl->GetAliRun() == 0x0) {
 121       rl->LoadgAlice();
 122       gAlice = rl->GetAliRun(); // should be removed
 123       rl->LoadHeader();
 124     }
 125   }
 126   fgPhosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
 127   if ( !fgPhosLoader )
 128     Error("AliPHOSGetter", "Could not find PHOSLoader") ;
 129   else
 130     fgPhosLoader->SetTitle(version);
 131
 132   // initialize data members
 133   SetDebug(0) ;
 134   fBTE = 0 ;
 135   fPrimaries = 0 ;
 136   fLoadingStatus = "" ;
 137
 138   fESDFileName = rl->GetFileName()  ; // this should be the galice.root file
 139   fESDFileName.ReplaceAll("galice.root", "AliESDs.root") ;
 140   fESDFile = 0 ;
 141   fESD = 0 ;
 142   fESDTree = 0 ;
 143   fRawDigits = kFALSE ;
 144
 145 }
 146
 147 AliPHOSGetter::AliPHOSGetter(const AliPHOSGetter & obj) :
 148   TObject(obj),
 149   fBTE(0),
 150   fLoadingStatus(),
 151   fNPrimaries(0),
 152   fPrimaries(0),
 153   fESDFile(0),
 154   fESDFileName(),
 155   fESD(0),
 156   fESDTree(0),
 157   fRawDigits(kFALSE),
 158   fcdb(0)
 159 {
 160   // cpy ctor requested by Coding Convention
 161   Fatal("cpy ctor", "not implemented") ;
 162 }
 163
 164 //____________________________________________________________________________
 165 AliPHOSGetter::AliPHOSGetter(Int_t /*i*/) :
 166   fBTE(0),
 167   fLoadingStatus(),
 168   fNPrimaries(0),
 169   fPrimaries(0),
 170   fESDFile(0),
 171   fESDFileName(),
 172   fESD(0),
 173   fESDTree(0),
 174   fRawDigits(kFALSE),
 175   fcdb(0)
 176 {
 177   // special constructor for onflight
 178 }
 179
 180
 181 //____________________________________________________________________________
 182 AliPHOSGetter::~AliPHOSGetter()
 183 {
 184   // dtor
 185   if(fgPhosLoader){
 186     delete fgPhosLoader ;
 187     fgPhosLoader = 0 ;
 188   }
 189   if(fBTE){
 190     delete fBTE ;
 191     fBTE = 0 ;
 192   }
 193   if(fPrimaries){
 194     fPrimaries->Delete() ;
 195     delete fPrimaries ;
 196   }
 197   if (fESD)
 198     delete fESD ;
 199   if (fESDTree)
 200     delete fESDTree ;
 201
 202   fgObjGetter = 0;
 203 }
 204
 205 //____________________________________________________________________________
 206 void AliPHOSGetter::Reset()
 207 {
 208   // resets things in case the getter is called consecutively with different files
 209   // the PHOS Loader is already deleted by the Run Loader
 210
 211   if (fPrimaries) {
 212     fPrimaries->Delete() ;
 213     delete fPrimaries ;
 214   }
 215   fgPhosLoader = 0;
 216   fgObjGetter = 0;
 217 }
 218
 219 //____________________________________________________________________________
 220 AliPHOSClusterizer * AliPHOSGetter::Clusterizer()
 221 {
 222   // Returns pointer to the Clusterizer task
 223   AliPHOSClusterizer * rv ;
 224   rv =  dynamic_cast<AliPHOSClusterizer *>(PhosLoader()->Reconstructioner()) ;
 225   if (!rv) {
 226     Event(0, "R") ;
 227     rv =  dynamic_cast<AliPHOSClusterizer*>(PhosLoader()->Reconstructioner()) ;
 228   }
 229   return rv ;
 230 }
 231
 232 //____________________________________________________________________________
 233 TObjArray * AliPHOSGetter::CpvRecPoints() const
 234 {
 235   // asks the Loader to return the CPV RecPoints container
 236
 237   TObjArray * rv = 0 ;
 238
 239   rv = PhosLoader()->CpvRecPoints() ;
 240   if (!rv) {
 241     PhosLoader()->MakeRecPointsArray() ;
 242     rv = PhosLoader()->CpvRecPoints() ;
 243   }
 244   return rv ;
 245 }
 246
 247 //____________________________________________________________________________
 248 TClonesArray * AliPHOSGetter::Digits() const
 249 {
 250   // asks the Loader to return the Digits container
 251
 252   TClonesArray * rv = 0 ;
 253   rv = PhosLoader()->Digits() ;
 254
 255   if( !rv ) {
 256     PhosLoader()->MakeDigitsArray() ;
 257     rv = PhosLoader()->Digits() ;
 258   }
 259   return rv ;
 260 }
 261
 262 //____________________________________________________________________________
 263 AliPHOSDigitizer * AliPHOSGetter::Digitizer()
 264 {
 265   // Returns pointer to the Digitizer task
 266   AliPHOSDigitizer * rv ;
 267   rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
 268   if (!rv) {
 269     Event(0, "D") ;
 270     rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
 271   }
 272   return rv ;
 273 }
 274
 275
 276 //____________________________________________________________________________
 277 TObjArray * AliPHOSGetter::EmcRecPoints() const
 278 {
 279   // asks the Loader to return the EMC RecPoints container
 280
 281   TObjArray * rv = 0 ;
 282
 283   rv = PhosLoader()->EmcRecPoints() ;
 284   if (!rv) {
 285     PhosLoader()->MakeRecPointsArray() ;
 286     rv = PhosLoader()->EmcRecPoints() ;
 287   }
 288   return rv ;
 289 }
 290
 291 //____________________________________________________________________________
 292 TClonesArray * AliPHOSGetter::TrackSegments() const
 293 {
 294   // asks the Loader to return the TrackSegments container
 295
 296   TClonesArray * rv = 0 ;
 297
 298   rv = PhosLoader()->TrackSegments() ;
 299   if (!rv) {
 300     PhosLoader()->MakeTrackSegmentsArray() ;
 301     rv = PhosLoader()->TrackSegments() ;
 302   }
 303   return rv ;
 304 }
 305
 306 //____________________________________________________________________________
 307 AliPHOSTrackSegmentMaker * AliPHOSGetter::TrackSegmentMaker()
 308 {
 309   // Returns pointer to the TrackSegmentMaker task
 310   AliPHOSTrackSegmentMaker * rv ;
 311   rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
 312   if (!rv) {
 313     Event(0, "T") ;
 314     rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
 315   }
 316   return rv ;
 317 }
 318
 319 //____________________________________________________________________________
 320 TClonesArray * AliPHOSGetter::RecParticles() const
 321 {
 322   // asks the Loader to return the TrackSegments container
 323
 324   TClonesArray * rv = 0 ;
 325
 326   rv = PhosLoader()->RecParticles() ;
 327   if (!rv) {
 328     PhosLoader()->MakeRecParticlesArray() ;
 329     rv = PhosLoader()->RecParticles() ;
 330   }
 331   return rv ;
 332 }
 333 //____________________________________________________________________________
 334 void AliPHOSGetter::Event(Int_t event, const char* opt)
 335 {
 336   // Reads the content of all Tree's S, D and R
 337
 338 //   if ( event >= MaxEvent() ) {
 339 //     Error("Event", "%d not found in TreeE !", event) ;
 340 //     return ;
 341 //   }
 342
 343   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 344
 345 //   // checks if we are dealing with test-beam data
 346 //   TBranch * btb = rl->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
 347 //   if(btb){
 348 //     if(!fBTE)
 349 //       fBTE = new AliPHOSBeamTestEvent() ;
 350 //     btb->SetAddress(&fBTE) ;
 351 //     btb->GetEntry(event) ;
 352 //   }
 353 //   else{
 354 //     if(fBTE){
 355 //       delete fBTE ;
 356 //       fBTE = 0 ;
 357 //     }
 358 //   }
 359
 360   // Loads the type of object(s) requested
 361
 362   rl->GetEvent(event) ;
 363
 364   if(strstr(opt,"X") || (strcmp(opt,"")==0)){
 365     ReadPrimaries() ;
 366   }
 367
 368   if(strstr(opt,"H")  ){
 369     ReadTreeH();
 370   }
 371
 372   if(strstr(opt,"S")  ){
 373     ReadTreeS() ;
 374   }
 375
 376   if(strstr(opt,"D") ){
 377     ReadTreeD() ;
 378   }
 379
 380   if(strstr(opt,"R") ){
 381     ReadTreeR() ;
 382   }
 383
 384   if( strstr(opt,"T") ){
 385     ReadTreeT() ;
 386   }
 387
 388   if( strstr(opt,"P") ){
 389     ReadTreeP() ;
 390   }
 391
 392   if( strstr(opt,"E") ){
 393     ReadTreeE(event) ;
 394   }
 395
 396 }
 397
 398
 399 //____________________________________________________________________________
 400 void AliPHOSGetter::Event(AliRawReader *rawReader, const char* opt, Bool_t isOldRCUFormat)
 401 {
 402   // Reads the raw event from rawReader
 403   // isOldRCUFormat defines whenever to assume
 404   // the old RCU format or not
 405
 406   if( strstr(opt,"W")  ){
 407     rawReader->NextEvent();
 408     ReadRaw(rawReader,isOldRCUFormat) ;
 409   }
 410
 411 }
 412
 413
 414 //____________________________________________________________________________
 415 Int_t AliPHOSGetter::EventNumber() const
 416   {
 417   // return the current event number
 418   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 419   return static_cast<Int_t>(rl->GetEventNumber()) ;
 420 }
 421
 422 //____________________________________________________________________________
 423   TClonesArray * AliPHOSGetter::Hits() const
 424 {
 425   // asks the loader to return  the Hits container
 426
 427   TClonesArray * rv = 0 ;
 428
 429   rv = PhosLoader()->Hits() ;
 430   if ( !rv ) {
 431     PhosLoader()->LoadHits("read");
 432     rv = PhosLoader()->Hits() ;
 433   }
 434   return rv ;
 435 }
 436
 437 //____________________________________________________________________________
 438 AliPHOSGetter * AliPHOSGetter::Instance(const char* alirunFileName, const char* version, Option_t * openingOption)
 439 {
 440   // Creates and returns the pointer of the unique instance
 441   // Must be called only when the environment has changed
 442
 443   if(!fgObjGetter){ // first time the getter is called
 444     fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 445   }
 446   else { // the getter has been called previously
 447     AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
 448     if ( rl->GetFileName() == alirunFileName ) {// the alirunFile has the same name
 449       // check if the file is already open
 450       TFile * galiceFile = dynamic_cast<TFile *>(gROOT->FindObject(rl->GetFileName()) ) ;
 451
 452       if ( !galiceFile )
 453         fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 454
 455       else {  // the file is already open check the version name
 456         TString currentVersionName = rl->GetEventFolder()->GetName() ;
 457         TString newVersionName(version) ;
 458         if (currentVersionName == newVersionName)
 459           if(fgDebug)
 460             ::Warning( "Instance", "Files with version %s already open", currentVersionName.Data() ) ;
 461         else {
 462           fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 463         }
 464       }
 465     }
 466     else {
 467       AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle()) ;
 468       if ( strstr(version, AliConfig::GetDefaultEventFolderName()) ) // false in case of merging
 469         delete rl ;
 470       fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
 471     }
 472   }
 473   if (!fgObjGetter)
 474     ::Error("AliPHOSGetter::Instance", "Failed to create the PHOS Getter object") ;
 475   else
 476     if (fgDebug)
 477       Print() ;
 478
 479   return fgObjGetter ;
 480 }
 481
 482 //____________________________________________________________________________
 483 AliPHOSGetter *  AliPHOSGetter::Instance()
 484 {
 485   // Returns the pointer of the unique instance already defined
 486
 487   if(!fgObjGetter && fgDebug)
 488      ::Warning("AliPHOSGetter::Instance", "Getter not initialized") ;
 489
 490    return fgObjGetter ;
 491
 492 }
 493
 494 //____________________________________________________________________________
 495 Int_t AliPHOSGetter::MaxEvent() const
 496 {
 497   // returns the number of events in the run (from TE)
 498
 499   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 500   return static_cast<Int_t>(rl->GetNumberOfEvents()) ;
 501 }
 502
 503 //____________________________________________________________________________
 504 TParticle * AliPHOSGetter::Primary(Int_t index) const
 505 {
 506   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 507   return rl->Stack()->Particle(index) ;
 508 }
 509
 510 //____________________________________________________________________________
 511 AliPHOS * AliPHOSGetter:: PHOS() const
 512 {
 513   // returns the PHOS object
 514   AliPHOSLoader *    loader = 0;
 515   static AliPHOSLoader * oldloader = 0;
 516   static AliPHOS * phos = 0;
 517
 518   loader = PhosLoader();
 519
 520   if ( loader != oldloader) {
 521     phos = dynamic_cast<AliPHOS*>(loader->GetModulesFolder()->FindObject("PHOS")) ;
 522     oldloader = loader;
 523   }
 524   if (!phos)
 525     if (fgDebug)
 526       Warning("PHOS", "PHOS module not found in module folders: %s", PhosLoader()->GetModulesFolder()->GetName() ) ;
 527   return phos ;
 528 }
 529
 530
 531
 532 //____________________________________________________________________________
 533 AliPHOSPID * AliPHOSGetter::PID()
 534 {
 535   // Returns pointer to the PID task
 536   AliPHOSPID * rv ;
 537   rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
 538   if (!rv) {
 539     Event(0, "P") ;
 540     rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
 541   }
 542   return rv ;
 543 }
 544
 545 //____________________________________________________________________________
 546 AliPHOSGeometry * AliPHOSGetter::PHOSGeometry() const
 547 {
 548   // Returns PHOS geometry
 549
 550   AliPHOSGeometry * rv = 0 ;
 551   if (PHOS() )
 552     rv =  PHOS()->GetGeometry() ;
 553   else {
 554     rv = AliPHOSGeometry::GetInstance();
 555     if (!rv) {
 556       AliError("Could not find PHOS geometry! Loading the default one !");
 557       rv = AliPHOSGeometry::GetInstance("IHEP","");
 558     }
 559   }
 560   return rv ;
 561 }
 562
 563 //____________________________________________________________________________
 564 TClonesArray * AliPHOSGetter::Primaries()
 565 {
 566   // creates the Primaries container if needed
 567   if ( !fPrimaries ) {
 568     if (fgDebug)
 569       Info("Primaries", "Creating a new TClonesArray for primaries") ;
 570     fPrimaries = new TClonesArray("TParticle", 1000) ;
 571   }
 572   return fPrimaries ;
 573 }
 574
 575 //____________________________________________________________________________
 576 void  AliPHOSGetter::Print()
 577 {
 578   // Print usefull information about the getter
 579
 580   AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
 581   ::Info( "Print", "gAlice file is %s -- version name is %s", (rl->GetFileName()).Data(), rl->GetEventFolder()->GetName() ) ;
 582 }
 583
 584 //____________________________________________________________________________
 585 void AliPHOSGetter::ReadPrimaries()
 586 {
 587   // Read Primaries from Kinematics.root
 588
 589   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
 590
 591   // gets kine tree from the root file (Kinematics.root)
 592   if ( ! rl->TreeK() ) { // load treeK the first time
 593     rl->LoadKinematics() ;
 594   }
 595
 596   fNPrimaries = (rl->GetHeader())->GetNtrack();
 597   if (fgDebug)
 598     Info( "ReadTreeK", "Found %d particles in event # %d", fNPrimaries, EventNumber() ) ;
 599
 600
 601   // first time creates the container
 602   if ( Primaries() )
 603     fPrimaries->Clear() ;
 604
 605   Int_t index = 0 ;
 606   for (index = 0 ; index < fNPrimaries; index++) {
 607     new ((*fPrimaries)[index]) TParticle(*(Primary(index)));
 608   }
 609 }
 610
 611 //____________________________________________________________________________
 612 Bool_t AliPHOSGetter::OpenESDFile()
 613 {
 614   //Open the ESD file
 615   Bool_t rv = kTRUE ;
 616   if (!fESDFile) {
 617     fESDFile = TFile::Open(fESDFileName) ;
 618     if (!fESDFile )
 619       return kFALSE ;
 620   }
 621   else if (fESDFile->IsOpen()) {
 622     fESDFile->Close() ;
 623     fESDFile = TFile::Open(fESDFileName) ;
 624   }
 625   if (!fESDFile->IsOpen())
 626     rv = kFALSE ;
 627   return rv ;
 628 }
 629
 630 //____________________________________________________________________________
 631 void AliPHOSGetter::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Double_t & energy, Double_t & time) const
 632 {
 633   // Fits the raw signal time distribution
 634
 635   const Int_t kNoiseThreshold = 0 ;
 636   Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
 637   Double_t signal = 0., signalmax = 0. ;
 638   time   = 0. ;
 639   energy = 0. ;
 640
 641   if (lowGainFlag) {
 642     timezero1 = timezero2 = signalmax = timemax = 0. ;
 643     signalF->FixParameter(0, PHOS()->GetRawFormatLowCharge()) ;
 644     signalF->FixParameter(1, PHOS()->GetRawFormatLowGain()) ;
 645     Int_t index ;
 646     for (index = 0; index < PHOS()->GetRawFormatTimeBins(); index++) {
 647       gLowGain->GetPoint(index, time, signal) ;
 648       if (signal > kNoiseThreshold && timezero1 == 0.)
 649         timezero1 = time ;
 650       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 651         timezero2 = time ;
 652       if (signal > signalmax) {
 653         signalmax = signal ;
 654         timemax   = time ;
 655       }
 656     }
 657     signalmax /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatLowCharge(),
 658                                                 PHOS()->GetRawFormatLowGain()) ;
 659     if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->GetRawFormatTimeMax() * 0.4 ) { // else its noise
 660       signalF->SetParameter(2, signalmax) ;
 661       signalF->SetParameter(3, timezero1) ;
 662       gLowGain->Fit(signalF, "QRO", "", 0., timezero2); //, "QRON") ;
 663       energy = signalF->GetParameter(2) ;
 664       time   = signalF->GetMaximumX() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
 665     }
 666   } else {
 667     timezero1 = timezero2 = signalmax = timemax = 0. ;
 668     signalF->FixParameter(0, PHOS()->GetRawFormatHighCharge()) ;
 669     signalF->FixParameter(1, PHOS()->GetRawFormatHighGain()) ;
 670     Int_t index ;
 671     for (index = 0; index < PHOS()->GetRawFormatTimeBins(); index++) {
 672       gHighGain->GetPoint(index, time, signal) ;
 673       if (signal > kNoiseThreshold && timezero1 == 0.)
 674         timezero1 = time ;
 675       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
 676         timezero2 = time ;
 677       if (signal > signalmax) {
 678         signalmax = signal ;
 679         timemax   = time ;
 680       }
 681     }
 682     signalmax /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatHighCharge(),
 683                                                 PHOS()->GetRawFormatHighGain()) ;;
 684     if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->GetRawFormatTimeMax() * 0.4 ) { // else its noise
 685       signalF->SetParameter(2, signalmax) ;
 686       signalF->SetParameter(3, timezero1) ;
 687       gHighGain->Fit(signalF, "QRO", "", 0., timezero2) ;
 688       energy = signalF->GetParameter(2) ;
 689       time   = signalF->GetMaximumX() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
 690     }
 691   }
 692   if (time == 0) energy = 0 ;
 693 }
 694
 695 //____________________________________________________________________________
 696 Int_t AliPHOSGetter::CalibrateRaw(Double_t energy, Int_t *relId)
 697 {
 698   // Convert energy into digitized amplitude for a cell relId
 699   // It is a user responsilibity to open CDB and set
 700   // AliPHOSCalibData object by the following operators:
 701   //
 702   // AliCDBLocal *loc = new AliCDBLocal("deCalibDB");
 703   // AliPHOSCalibData* clb = (AliPHOSCalibData*)AliCDBStorage::Instance()
 704   //    ->Get(path_to_calibdata,run_number);
 705   // AliPHOSGetter* gime = AliPHOSGetter::Instance("galice.root");
 706   // gime->SetCalibData(clb);
 707
 708   if (CalibData() == 0)
 709     Warning("CalibrateRaw","Calibration DB is not initiated!");
 710
 711   Int_t   module = relId[0];
 712   Int_t   column = relId[3];
 713   Int_t   row    = relId[2];
 714
 715   Float_t gainFactor = 0.0015; // width of one Emc ADC channel in GeV
 716   Float_t pedestal   = 0.005;  // Emc pedestals
 717
 718   if(CalibData()) {
 719     gainFactor = CalibData()->GetADCchannelEmc (module,column,row);
 720     pedestal   = CalibData()->GetADCpedestalEmc(module,column,row);
 721   }
 722
 723   Int_t   amp = static_cast<Int_t>( (energy - pedestal) / gainFactor + 0.5 ) ;
 724   return amp;
 725 }
 726 //____________________________________________________________________________
 727 Int_t AliPHOSGetter::ReadRaw(AliRawReader *rawReader,Bool_t isOldRCUFormat)
 728 {
 729   // reads the raw format data, converts it into digits format and store digits in Digits()
 730   // container.
 731   // isOldRCUFormat = kTRUE in case of the old RCU
 732   // format used in the raw data readout.
 733   // Reimplemented by Boris Polichtchouk (Jul 2006)
 734   // to make it working with the Jul-Aug 2006 beam test data.
 735
 736   AliCaloRawStream in(rawReader,"PHOS");
 737   in.SetOldRCUFormat(isOldRCUFormat);
 738
 739   TF1 * signalF = new TF1("signal", AliPHOS::RawResponseFunction, 0, PHOS()->GetRawFormatTimeMax(), 4);
 740   signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ;
 741
 742   Int_t relId[4], absId =0;
 743   Bool_t lowGainFlag = kFALSE ;
 744
 745   TClonesArray * digits = Digits() ;
 746   digits->Clear() ;
 747
 748   Int_t    iBin     = 0;
 749   Int_t    iDigit   = 0 ;
 750   Double_t energyHG = 0. ;
 751   Double_t energyLG = 0. ;
 752   Double_t time     = 0. ;
 753   Int_t    iOldDigit;
 754   Bool_t   seen;
 755   TH1F*    hLowGain  = 0;
 756   TH1F*    hHighGain = 0;
 757
 758   while ( in.Next() ) { // PHOS entries loop
 759
 760     lowGainFlag = in.IsLowGain();
 761
 762     // (Re)create histograms to store samples
 763
 764     if (lowGainFlag) {
 765       if(!hLowGain)
 766         hLowGain = new TH1F("hLowGain","Low gain",in.GetTimeLength(),0,in.GetTimeLength());
 767       else
 768         if(hLowGain->GetNbinsX() != in.GetTimeLength()) {
 769           delete hLowGain;
 770           hLowGain = new TH1F("hLowGain","Low gain",in.GetTimeLength(),0,in.GetTimeLength());
 771         }
 772     } else {
 773       if(!hHighGain)
 774         hHighGain = new TH1F("hHighGain","High gain",in.GetTimeLength(),0,in.GetTimeLength());
 775       else
 776         if(hHighGain->GetNbinsX() != in.GetTimeLength()) {
 777           delete hHighGain;
 778           hHighGain = new TH1F("hHighGain","High gain",in.GetTimeLength(),0,in.GetTimeLength());
 779         }
 780     }
 781
 782     // Fill histograms with samples
 783     if(lowGainFlag)
 784       hLowGain ->SetBinContent(in.GetTimeLength()-iBin-1,in.GetSignal());
 785     else {
 786       hHighGain->SetBinContent(in.GetTimeLength()-iBin-1,in.GetSignal());
 787     }
 788
 789     iBin++;
 790
 791     // Fit the full sample
 792     if(iBin==in.GetTimeLength()) {
 793       iBin=0;
 794
 795       // Temporarily we do not fit the sample graph, but
 796       // take the energy from the graph maximum, and the pedestal from the 0th point
 797       // 30 Aug 2006
 798
 799       //FitRaw(lowGainFlag, hLowGain, hHighGain, signalF, energy, time);
 800
 801       if(lowGainFlag) {
 802         energyLG  = hLowGain ->GetMaximum();     // "digit amplitude"
 803 //      energyLG -= hLowGain ->GetBinContent(0); // "pedestal subtraction"
 804         energyLG *= AliPHOS::GetRawFormatHighLowGainFactor(); // *16
 805         if(AliLog::GetGlobalDebugLevel()>3)
 806           AliDebug(4,Form("----Printing hLowGain: ----\n")) ; hLowGain ->Print("all");
 807         AliDebug(2,Form("AliPHOSGetter::ReadRaw: (mod,col,row)=(%d,%d,%d), low gain energy=%f\n\n",
 808                         in.GetModule(),in.GetColumn(),in.GetRow(),energyLG));
 809       }
 810
 811       else {
 812         energyHG  = hHighGain->GetMaximum();     // "digit amplitude"
 813 //      energyHG -= hHighGain->GetBinContent(0); // "pedestal subtraction"
 814         if(AliLog::GetGlobalDebugLevel()>3)
 815           AliDebug(4,Form("----Printing hHighGain: ----\n")); hHighGain->Print("all");
 816         AliDebug(2,Form("AliPHOSGetter::ReadRaw: (mod,col,row)=(%d,%d,%d), high gain energy=%f\n\n",
 817                         in.GetModule(),in.GetColumn(),in.GetRow(),energyHG));
 818       }
 819
 820       // Time is not evaluated for the moment (12.01.2007). To be implemented later.
 821       time = -1;
 822
 823       relId[0] = in.GetModule()+1;
 824       relId[1] = 0;
 825       relId[2] = in.GetRow();
 826       relId[3] = in.GetColumn();
 827       if(!PHOSGeometry()) AliFatal("Couldn't find PHOSGeometry!");
 828       PHOSGeometry()->RelToAbsNumbering(relId, absId);
 829
 830       // Add low gain digit only if the high gain digit does not exist in the digits array
 831       if(lowGainFlag) {
 832         seen = kFALSE;
 833         for (iOldDigit=iDigit; iOldDigit=0; iOldDigit--) {
 834           if (dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit))->GetId() == absId) {
 835             seen = kTRUE;
 836             break;
 837           }
 838         }
 839         if (!seen)
 840           new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyLG,time);
 841         energyLG = 0. ;
 842       }
 843       // Add high gain digit only if the low gain digit does not exist in the digits array
 844       // and it is not saturated
 845       else {
 846         seen = kFALSE;
 847         for (iOldDigit=iDigit; iOldDigit=0; iOldDigit--) {
 848           if (dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit))->GetId() == absId) {
 849             seen = kTRUE;
 850             break;
 851           }
 852         }
 853         if (!seen && energyHG < 1023)
 854           new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyHG,time);
 855         energyHG = 0. ;
 856       }
 857       iDigit++;
 858     }
 859   }
 860
 861   // PHOS entries loop
 862
 863   digits->Sort() ;
 864
 865   //!!!!for debug!!!
 866   Int_t modMax=-111;
 867   Int_t colMax=-111;
 868   Int_t rowMax=-111;
 869   Float_t eMax=-333;
 870   //!!!for debug!!!
 871
 872   for(iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
 873     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
 874     if(digit->GetEnergy()>eMax) {
 875       PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
 876       eMax=digit->GetEnergy();
 877       modMax=relId[0];
 878       rowMax=relId[2];
 879       colMax=relId[3];
 880     }
 881   }
 882
 883   AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
 884                   modMax,colMax,rowMax,eMax));
 885
 886   delete signalF ;
 887   delete hHighGain;
 888   delete hLowGain;
 889
 890   return digits->GetEntriesFast() ;
 891 }
 892
 893 //____________________________________________________________________________
 894 Int_t AliPHOSGetter::ReadTreeD() const
 895 {
 896   // Read the Digits
 897
 898   PhosLoader()->CleanDigits() ;
 899   PhosLoader()->LoadDigits("UPDATE") ;
 900   PhosLoader()->LoadDigitizer("UPDATE") ;
 901   return Digits()->GetEntries() ;
 902 }
 903
 904 //____________________________________________________________________________
 905 Int_t AliPHOSGetter::ReadTreeH() const
 906 {
 907   // Read the Hits
 908   PhosLoader()->CleanHits() ;
 909   // gets TreeH from the root file (PHOS.Hit.root)
 910   //if ( !IsLoaded("H") ) {
 911     PhosLoader()->LoadHits("UPDATE") ;
 912   //  SetLoaded("H") ;
 913   //}
 914   return Hits()->GetEntries() ;
 915 }
 916
 917 //____________________________________________________________________________
 918 Int_t AliPHOSGetter::ReadTreeR() const
 919 {
 920   // Read the RecPoints
 921
 922   PhosLoader()->CleanRecPoints() ;
 923   // gets TreeR from the root file (PHOS.RecPoints.root)
 924   //if ( !IsLoaded("R") ) {
 925     PhosLoader()->LoadRecPoints("UPDATE") ;
 926     PhosLoader()->LoadClusterizer("UPDATE") ;
 927     //  SetLoaded("R") ;
 928     //}
 929
 930   return EmcRecPoints()->GetEntries() ;
 931 }
 932
 933 //____________________________________________________________________________
 934 Int_t AliPHOSGetter::ReadTreeT() const
 935 {
 936   // Read the TrackSegments
 937
 938   PhosLoader()->CleanTracks() ;
 939   // gets TreeT from the root file (PHOS.TrackSegments.root)
 940   //if ( !IsLoaded("T") ) {
 941     PhosLoader()->LoadTracks("UPDATE") ;
 942     PhosLoader()->LoadTrackSegmentMaker("UPDATE") ;
 943     //    SetLoaded("T") ;
 944     //}
 945
 946   return TrackSegments()->GetEntries() ;
 947 }
 948 //____________________________________________________________________________
 949 Int_t AliPHOSGetter::ReadTreeP() const
 950 {
 951   // Read the RecParticles
 952
 953   PhosLoader()->CleanRecParticles() ;
 954
 955   // gets TreeT from the root file (PHOS.TrackSegments.root)
 956   //  if ( !IsLoaded("P") ) {
 957     PhosLoader()->LoadRecParticles("UPDATE") ;
 958     PhosLoader()->LoadPID("UPDATE") ;
 959     //  SetLoaded("P") ;
 960     //}
 961
 962   return RecParticles()->GetEntries() ;
 963 }
 964 //____________________________________________________________________________
 965 Int_t AliPHOSGetter::ReadTreeS() const
 966 {
 967   // Read the SDigits
 968
 969   PhosLoader()->CleanSDigits() ;
 970   // gets TreeS from the root file (PHOS.SDigits.root)
 971   //if ( !IsLoaded("S") ) {
 972     PhosLoader()->LoadSDigits("READ") ;
 973     PhosLoader()->LoadSDigitizer("READ") ;
 974     //  SetLoaded("S") ;
 975     //}
 976
 977   return SDigits()->GetEntries() ;
 978 }
 979
 980 //____________________________________________________________________________
 981 Int_t AliPHOSGetter::ReadTreeE(Int_t event)
 982 {
 983   // Read the ESD
 984
 985   // gets esdTree from the root file (AliESDs.root)
 986   if (!fESDFile)
 987     if ( !OpenESDFile() )
 988       return -1 ;
 989
 990   fESDTree = static_cast<TTree*>(fESDFile->Get("esdTree")) ;
 991   fESD = new AliESD;
 992    if (!fESDTree) {
 993
 994      Error("ReadTreeE", "no ESD tree found");
 995      return -1;
 996    }
 997    fESDTree->SetBranchAddress("ESD", &fESD);
 998    fESDTree->GetEvent(event);
 999
1000    return event ;
1001 }
1002
1003 //____________________________________________________________________________
1004 TClonesArray * AliPHOSGetter::SDigits() const
1005 {
1006   // asks the Loader to return the Digits container
1007
1008   TClonesArray * rv = 0 ;
1009
1010   rv = PhosLoader()->SDigits() ;
1011   if (!rv) {
1012     PhosLoader()->MakeSDigitsArray() ;
1013     rv = PhosLoader()->SDigits() ;
1014   }
1015   return rv ;
1016 }
1017
1018 //____________________________________________________________________________
1019 AliPHOSSDigitizer * AliPHOSGetter::SDigitizer()
1020 {
1021   // Returns pointer to the SDigitizer task
1022   AliPHOSSDigitizer * rv ;
1023   rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1024   if (!rv) {
1025     Event(0, "S") ;
1026     rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1027   }
1028   return rv ;
1029 }
1030
1031 //____________________________________________________________________________
1032 TParticle * AliPHOSGetter::Secondary(const TParticle* p, Int_t index) const
1033 {
1034   // Return first (index=1) or second (index=2) secondary particle of primary particle p
1035
1036   if(index <= 0)
1037     return 0 ;
1038   if(index > 2)
1039     return 0 ;
1040
1041   if(p) {
1042   Int_t daughterIndex = p->GetDaughter(index-1) ;
1043   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1044   return  rl->GetAliRun()->GetMCApp()->Particle(daughterIndex) ;
1045   }
1046   else
1047     return 0 ;
1048 }
1049
1050 //____________________________________________________________________________
1051 void AliPHOSGetter::Track(Int_t itrack)
1052 {
1053   // Read the first entry of PHOS branch in hit tree gAlice->TreeH()
1054
1055  AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1056
1057   if( !TreeH() ) // load treeH the first time
1058     rl->LoadHits() ;
1059
1060   // first time create the container
1061   TClonesArray * hits = Hits() ;
1062   if ( hits )
1063     hits->Clear() ;
1064
1065   TBranch * phosbranch = dynamic_cast<TBranch*>(TreeH()->GetBranch("PHOS")) ;
1066   phosbranch->SetAddress(&hits) ;
1067   phosbranch->GetEntry(itrack) ;
1068 }
1069
1070 //____________________________________________________________________________
1071 TTree * AliPHOSGetter::TreeD() const
1072 {
1073   // Returns pointer to the Digits Tree
1074   TTree * rv = 0 ;
1075   rv = PhosLoader()->TreeD() ;
1076   if ( !rv ) {
1077     PhosLoader()->MakeTree("D");
1078     rv = PhosLoader()->TreeD() ;
1079   }
1080
1081   return rv ;
1082 }
1083
1084 //____________________________________________________________________________
1085 TTree * AliPHOSGetter::TreeH() const
1086 {
1087   // Returns pointer to the Hits Tree
1088   TTree * rv = 0 ;
1089   rv = PhosLoader()->TreeH() ;
1090   if ( !rv ) {
1091     PhosLoader()->MakeTree("H");
1092     rv = PhosLoader()->TreeH() ;
1093   }
1094
1095   return rv ;
1096 }
1097
1098 //____________________________________________________________________________
1099 TTree * AliPHOSGetter::TreeR() const
1100 {
1101   // Returns pointer to the RecPoints Tree
1102   TTree * rv = 0 ;
1103   rv = PhosLoader()->TreeR() ;
1104   if ( !rv ) {
1105     PhosLoader()->MakeTree("R");
1106     rv = PhosLoader()->TreeR() ;
1107   }
1108
1109   return rv ;
1110 }
1111
1112 //____________________________________________________________________________
1113 TTree * AliPHOSGetter::TreeT() const
1114 {
1115   // Returns pointer to the TrackSegments Tree
1116   TTree * rv = 0 ;
1117   rv = PhosLoader()->TreeT() ;
1118   if ( !rv ) {
1119     PhosLoader()->MakeTree("T");
1120     rv = PhosLoader()->TreeT() ;
1121   }
1122
1123   return rv ;
1124 }
1125 //____________________________________________________________________________
1126 TTree * AliPHOSGetter::TreeP() const
1127 {
1128   // Returns pointer to the RecParticles  Tree
1129   TTree * rv = 0 ;
1130   rv = PhosLoader()->TreeP() ;
1131   if ( !rv ) {
1132     PhosLoader()->MakeTree("P");
1133     rv = PhosLoader()->TreeP() ;
1134   }
1135
1136   return rv ;
1137 }
1138
1139 //____________________________________________________________________________
1140 TTree * AliPHOSGetter::TreeS() const
1141 {
1142  // Returns pointer to the SDigits Tree
1143   TTree * rv = 0 ;
1144   rv = PhosLoader()->TreeS() ;
1145   if ( !rv ) {
1146     PhosLoader()->MakeTree("S");
1147     rv = PhosLoader()->TreeS() ;
1148   }
1149
1150   return rv ;
1151 }
1152
1153 //____________________________________________________________________________
1154 Bool_t AliPHOSGetter::VersionExists(TString & opt) const
1155 {
1156   // checks if the version with the present name already exists in the same directory
1157
1158   Bool_t rv = kFALSE ;
1159
1160   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1161   TString version( rl->GetEventFolder()->GetName() ) ;
1162
1163   opt.ToLower() ;
1164
1165   if ( opt == "sdigits") {
1166     // add the version name to the root file name
1167     TString fileName( PhosLoader()->GetSDigitsFileName() ) ;
1168     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1169       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1170     if ( !(gSystem->AccessPathName(fileName)) ) {
1171       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1172       rv = kTRUE ;
1173     }
1174     PhosLoader()->SetSDigitsFileName(fileName) ;
1175   }
1176
1177   if ( opt == "digits") {
1178     // add the version name to the root file name
1179     TString fileName( PhosLoader()->GetDigitsFileName() ) ;
1180     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1181       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1182     if ( !(gSystem->AccessPathName(fileName)) ) {
1183       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1184       rv = kTRUE ;
1185     }
1186   }
1187
1188   return rv ;
1189
1190 }
1191
1192 //____________________________________________________________________________
1193 UShort_t AliPHOSGetter::EventPattern(void) const
1194 {
1195   // Return the pattern (trigger bit register) of the beam-test event
1196   if(fBTE)
1197     return fBTE->GetPattern() ;
1198   else
1199     return 0 ;
1200 }
1201 //____________________________________________________________________________
1202 Float_t AliPHOSGetter::BeamEnergy(void) const
1203 {
1204   // Return the beam energy of the beam-test event
1205   if(fBTE)
1206     return fBTE->GetBeamEnergy() ;
1207   else
1208     return 0 ;
1209 }
1210 //____________________________________________________________________________
1211
1212 AliPHOSCalibData* AliPHOSGetter::CalibData()
1213 {
1214   // Check if the instance of AliPHOSCalibData exists, and return it
1215
1216   return fgCalibData;
1217 }