Extract ALTRO sample generation to a separate class AliPHOSPulseGenerator
[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 "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 }