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