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 = new AliPHOSPulseGenerator();
 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 = new AliPHOSPulseGenerator();
 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")) ; hLowGain ->Print("all");
 814         AliDebug(2,Form("AliPHOSGetter::ReadRaw: (mod,col,row)=(%d,%d,%d), low gain energy=%f\n\n",
 815                         in.GetModule(),in.GetColumn(),in.GetRow(),energyLG));
 816       }
 817
 818       else {
 819         energyHG  = hHighGain->GetMaximum();     // "digit amplitude"
 820 //      energyHG -= hHighGain->GetBinContent(0); // "pedestal subtraction"
 821         if(AliLog::GetGlobalDebugLevel()>3)
 822           AliDebug(4,Form("----Printing hHighGain: ----\n")); hHighGain->Print("all");
 823         AliDebug(2,Form("AliPHOSGetter::ReadRaw: (mod,col,row)=(%d,%d,%d), high gain energy=%f\n\n",
 824                         in.GetModule(),in.GetColumn(),in.GetRow(),energyHG));
 825       }
 826
 827       // Time is not evaluated for the moment (12.01.2007). To be implemented later.
 828       time = -1;
 829
 830       relId[0] = in.GetModule()+1;
 831       relId[1] = 0;
 832       relId[2] = in.GetRow();
 833       relId[3] = in.GetColumn();
 834       if(!PHOSGeometry()) AliFatal("Couldn't find PHOSGeometry!");
 835       PHOSGeometry()->RelToAbsNumbering(relId, absId);
 836
 837       // Add low gain digit only if the high gain digit does not exist in the digits array
 838       if(lowGainFlag) {
 839         seen = kFALSE;
 840         for (iOldDigit=iDigit; iOldDigit==0; iOldDigit--) {
 841           if (dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit))->GetId() == absId) {
 842             seen = kTRUE;
 843             break;
 844           }
 845         }
 846         if (!seen)
 847           new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyLG,time);
 848         energyLG = 0. ;
 849       }
 850       // Add high gain digit only if it is not saturated;
 851       // replace low gain digit by a high gain one
 852       else {
 853         if (energyHG >= 1023) continue;
 854         for (iOldDigit=iDigit; iOldDigit==0; iOldDigit--) {
 855           if (dynamic_cast<AliPHOSDigit*>(digits->At(iOldDigit))->GetId() == absId) {
 856             digits->RemoveAt(iOldDigit);
 857             new((*digits)[iDigit]) AliPHOSDigit(-1,absId,(Float_t)energyHG,time);
 858             break;
 859           }
 860         }
 861         energyHG = 0. ;
 862       }
 863       iDigit++;
 864     }
 865   }
 866
 867   // PHOS entries loop
 868
 869   digits->Sort() ;
 870
 871   //!!!!for debug!!!
 872   Int_t modMax=-111;
 873   Int_t colMax=-111;
 874   Int_t rowMax=-111;
 875   Float_t eMax=-333;
 876   //!!!for debug!!!
 877
 878   for(iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
 879     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
 880     if(digit->GetEnergy()>eMax) {
 881       PHOSGeometry()->AbsToRelNumbering(digit->GetId(),relId);
 882       eMax=digit->GetEnergy();
 883       modMax=relId[0];
 884       rowMax=relId[2];
 885       colMax=relId[3];
 886     }
 887   }
 888
 889   AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
 890                   modMax,colMax,rowMax,eMax));
 891
 892   delete signalF ;
 893   delete hHighGain;
 894   delete hLowGain;
 895
 896   return digits->GetEntriesFast() ;
 897 }
 898
 899 //____________________________________________________________________________
 900 Int_t AliPHOSGetter::ReadTreeD() const
 901 {
 902   // Read the Digits
 903
 904   PhosLoader()->CleanDigits() ;
 905   PhosLoader()->LoadDigits("UPDATE") ;
 906   PhosLoader()->LoadDigitizer("UPDATE") ;
 907   return Digits()->GetEntries() ;
 908 }
 909
 910 //____________________________________________________________________________
 911 Int_t AliPHOSGetter::ReadTreeH() const
 912 {
 913   // Read the Hits
 914   PhosLoader()->CleanHits() ;
 915   // gets TreeH from the root file (PHOS.Hit.root)
 916   //if ( !IsLoaded("H") ) {
 917     PhosLoader()->LoadHits("UPDATE") ;
 918   //  SetLoaded("H") ;
 919   //}
 920   return Hits()->GetEntries() ;
 921 }
 922
 923 //____________________________________________________________________________
 924 Int_t AliPHOSGetter::ReadTreeR() const
 925 {
 926   // Read the RecPoints
 927
 928   PhosLoader()->CleanRecPoints() ;
 929   // gets TreeR from the root file (PHOS.RecPoints.root)
 930   //if ( !IsLoaded("R") ) {
 931     PhosLoader()->LoadRecPoints("UPDATE") ;
 932     PhosLoader()->LoadClusterizer("UPDATE") ;
 933     //  SetLoaded("R") ;
 934     //}
 935
 936   return EmcRecPoints()->GetEntries() ;
 937 }
 938
 939 //____________________________________________________________________________
 940 Int_t AliPHOSGetter::ReadTreeT() const
 941 {
 942   // Read the TrackSegments
 943
 944   PhosLoader()->CleanTracks() ;
 945   // gets TreeT from the root file (PHOS.TrackSegments.root)
 946   //if ( !IsLoaded("T") ) {
 947     PhosLoader()->LoadTracks("UPDATE") ;
 948     PhosLoader()->LoadTrackSegmentMaker("UPDATE") ;
 949     //    SetLoaded("T") ;
 950     //}
 951
 952   return TrackSegments()->GetEntries() ;
 953 }
 954 //____________________________________________________________________________
 955 Int_t AliPHOSGetter::ReadTreeP() const
 956 {
 957   // Read the RecParticles
 958
 959   PhosLoader()->CleanRecParticles() ;
 960
 961   // gets TreeT from the root file (PHOS.TrackSegments.root)
 962   //  if ( !IsLoaded("P") ) {
 963     PhosLoader()->LoadRecParticles("UPDATE") ;
 964     PhosLoader()->LoadPID("UPDATE") ;
 965     //  SetLoaded("P") ;
 966     //}
 967
 968   return RecParticles()->GetEntries() ;
 969 }
 970 //____________________________________________________________________________
 971 Int_t AliPHOSGetter::ReadTreeS() const
 972 {
 973   // Read the SDigits
 974
 975   PhosLoader()->CleanSDigits() ;
 976   // gets TreeS from the root file (PHOS.SDigits.root)
 977   //if ( !IsLoaded("S") ) {
 978     PhosLoader()->LoadSDigits("READ") ;
 979     PhosLoader()->LoadSDigitizer("READ") ;
 980     //  SetLoaded("S") ;
 981     //}
 982
 983   return SDigits()->GetEntries() ;
 984 }
 985
 986 //____________________________________________________________________________
 987 Int_t AliPHOSGetter::ReadTreeE(Int_t event)
 988 {
 989   // Read the ESD
 990
 991   // gets esdTree from the root file (AliESDs.root)
 992   if (!fESDFile)
 993     if ( !OpenESDFile() )
 994       return -1 ;
 995
 996   fESDTree = static_cast<TTree*>(fESDFile->Get("esdTree")) ;
 997   fESD = new AliESD;
 998    if (!fESDTree) {
 999
1000      Error("ReadTreeE", "no ESD tree found");
1001      return -1;
1002    }
1003    fESDTree->SetBranchAddress("ESD", &fESD);
1004    fESDTree->GetEvent(event);
1005
1006    return event ;
1007 }
1008
1009 //____________________________________________________________________________
1010 TClonesArray * AliPHOSGetter::SDigits() const
1011 {
1012   // asks the Loader to return the Digits container
1013
1014   TClonesArray * rv = 0 ;
1015
1016   rv = PhosLoader()->SDigits() ;
1017   if (!rv) {
1018     PhosLoader()->MakeSDigitsArray() ;
1019     rv = PhosLoader()->SDigits() ;
1020   }
1021   return rv ;
1022 }
1023
1024 //____________________________________________________________________________
1025 AliPHOSSDigitizer * AliPHOSGetter::SDigitizer()
1026 {
1027   // Returns pointer to the SDigitizer task
1028   AliPHOSSDigitizer * rv ;
1029   rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1030   if (!rv) {
1031     Event(0, "S") ;
1032     rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
1033   }
1034   return rv ;
1035 }
1036
1037 //____________________________________________________________________________
1038 TParticle * AliPHOSGetter::Secondary(const TParticle* p, Int_t index) const
1039 {
1040   // Return first (index=1) or second (index=2) secondary particle of primary particle p
1041
1042   if(index <= 0)
1043     return 0 ;
1044   if(index > 2)
1045     return 0 ;
1046
1047   if(p) {
1048   Int_t daughterIndex = p->GetDaughter(index-1) ;
1049   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1050   return  rl->GetAliRun()->GetMCApp()->Particle(daughterIndex) ;
1051   }
1052   else
1053     return 0 ;
1054 }
1055
1056 //____________________________________________________________________________
1057 void AliPHOSGetter::Track(Int_t itrack)
1058 {
1059   // Read the first entry of PHOS branch in hit tree gAlice->TreeH()
1060
1061  AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1062
1063   if( !TreeH() ) // load treeH the first time
1064     rl->LoadHits() ;
1065
1066   // first time create the container
1067   TClonesArray * hits = Hits() ;
1068   if ( hits )
1069     hits->Clear() ;
1070
1071   TBranch * phosbranch = dynamic_cast<TBranch*>(TreeH()->GetBranch("PHOS")) ;
1072   phosbranch->SetAddress(&hits) ;
1073   phosbranch->GetEntry(itrack) ;
1074 }
1075
1076 //____________________________________________________________________________
1077 TTree * AliPHOSGetter::TreeD() const
1078 {
1079   // Returns pointer to the Digits Tree
1080   TTree * rv = 0 ;
1081   rv = PhosLoader()->TreeD() ;
1082   if ( !rv ) {
1083     PhosLoader()->MakeTree("D");
1084     rv = PhosLoader()->TreeD() ;
1085   }
1086
1087   return rv ;
1088 }
1089
1090 //____________________________________________________________________________
1091 TTree * AliPHOSGetter::TreeH() const
1092 {
1093   // Returns pointer to the Hits Tree
1094   TTree * rv = 0 ;
1095   rv = PhosLoader()->TreeH() ;
1096   if ( !rv ) {
1097     PhosLoader()->MakeTree("H");
1098     rv = PhosLoader()->TreeH() ;
1099   }
1100
1101   return rv ;
1102 }
1103
1104 //____________________________________________________________________________
1105 TTree * AliPHOSGetter::TreeR() const
1106 {
1107   // Returns pointer to the RecPoints Tree
1108   TTree * rv = 0 ;
1109   rv = PhosLoader()->TreeR() ;
1110   if ( !rv ) {
1111     PhosLoader()->MakeTree("R");
1112     rv = PhosLoader()->TreeR() ;
1113   }
1114
1115   return rv ;
1116 }
1117
1118 //____________________________________________________________________________
1119 TTree * AliPHOSGetter::TreeT() const
1120 {
1121   // Returns pointer to the TrackSegments Tree
1122   TTree * rv = 0 ;
1123   rv = PhosLoader()->TreeT() ;
1124   if ( !rv ) {
1125     PhosLoader()->MakeTree("T");
1126     rv = PhosLoader()->TreeT() ;
1127   }
1128
1129   return rv ;
1130 }
1131 //____________________________________________________________________________
1132 TTree * AliPHOSGetter::TreeP() const
1133 {
1134   // Returns pointer to the RecParticles  Tree
1135   TTree * rv = 0 ;
1136   rv = PhosLoader()->TreeP() ;
1137   if ( !rv ) {
1138     PhosLoader()->MakeTree("P");
1139     rv = PhosLoader()->TreeP() ;
1140   }
1141
1142   return rv ;
1143 }
1144
1145 //____________________________________________________________________________
1146 TTree * AliPHOSGetter::TreeS() const
1147 {
1148  // Returns pointer to the SDigits Tree
1149   TTree * rv = 0 ;
1150   rv = PhosLoader()->TreeS() ;
1151   if ( !rv ) {
1152     PhosLoader()->MakeTree("S");
1153     rv = PhosLoader()->TreeS() ;
1154   }
1155
1156   return rv ;
1157 }
1158
1159 //____________________________________________________________________________
1160 Bool_t AliPHOSGetter::VersionExists(TString & opt) const
1161 {
1162   // checks if the version with the present name already exists in the same directory
1163
1164   Bool_t rv = kFALSE ;
1165
1166   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
1167   TString version( rl->GetEventFolder()->GetName() ) ;
1168
1169   opt.ToLower() ;
1170
1171   if ( opt == "sdigits") {
1172     // add the version name to the root file name
1173     TString fileName( PhosLoader()->GetSDigitsFileName() ) ;
1174     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1175       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1176     if ( !(gSystem->AccessPathName(fileName)) ) {
1177       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1178       rv = kTRUE ;
1179     }
1180     PhosLoader()->SetSDigitsFileName(fileName) ;
1181   }
1182
1183   if ( opt == "digits") {
1184     // add the version name to the root file name
1185     TString fileName( PhosLoader()->GetDigitsFileName() ) ;
1186     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name
1187       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1188     if ( !(gSystem->AccessPathName(fileName)) ) {
1189       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1190       rv = kTRUE ;
1191     }
1192   }
1193
1194   return rv ;
1195
1196 }
1197
1198 //____________________________________________________________________________
1199 UShort_t AliPHOSGetter::EventPattern(void) const
1200 {
1201   // Return the pattern (trigger bit register) of the beam-test event
1202   if(fBTE)
1203     return fBTE->GetPattern() ;
1204   else
1205     return 0 ;
1206 }
1207 //____________________________________________________________________________
1208 Float_t AliPHOSGetter::BeamEnergy(void) const
1209 {
1210   // Return the beam energy of the beam-test event
1211   if(fBTE)
1212     return fBTE->GetBeamEnergy() ;
1213   else
1214     return 0 ;
1215 }
1216 //____________________________________________________________________________
1217
1218 AliPHOSCalibData* AliPHOSGetter::CalibData()
1219 {
1220   // Check if the instance of AliPHOSCalibData exists, and return it
1221
1222   return fgCalibData;
1223 }