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