Restoring raw data fit from version of 29-Aug-2004
[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 <TFrame.h>
51
52 // --- Standard library ---
53
54 // --- AliRoot header files ---
55 #include "AliESD.h"
56 #include "AliHeader.h"  
57 #include "AliMC.h"
58 #include "AliPHOS.h"
59 #include "AliPHOSBeamTestEvent.h"
60 #include "AliPHOSGetter.h"
61 #include "AliPHOSLoader.h"
62 #include "AliRunLoader.h"
63 #include "AliStack.h"  
64 #include "AliPHOSRawStream.h"
65 #include "AliRawReaderFile.h"
66 #include "AliLog.h"
67
68 ClassImp(AliPHOSGetter)
69   
70 AliPHOSGetter * AliPHOSGetter::fgObjGetter = 0 ; 
71 AliPHOSLoader * AliPHOSGetter::fgPhosLoader = 0;
72 Int_t AliPHOSGetter::fgDebug = 0;
73
74 //  TFile * AliPHOSGetter::fgFile = 0 ; 
75
76 //____________________________________________________________________________ 
77 AliPHOSGetter::AliPHOSGetter(const char* headerFile, const char* version, Option_t * openingOption)
78 {
79   // ctor only called by Instance()
80
81   AliRunLoader* rl = AliRunLoader::GetRunLoader(version) ; 
82   if (!rl) {
83     rl = AliRunLoader::Open(headerFile, version, openingOption);
84     if (!rl) {
85       Fatal("AliPHOSGetter", "Could not find the Run Loader for %s - %s",headerFile, version) ; 
86       return ;
87     } 
88     if (rl->GetAliRun() == 0x0) {
89       rl->LoadgAlice();
90       gAlice = rl->GetAliRun(); // should be removed
91       rl->LoadHeader();
92     }
93   }
94   fgPhosLoader = dynamic_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
95   if ( !fgPhosLoader ) 
96     Error("AliPHOSGetter", "Could not find PHOSLoader") ; 
97   else 
98     fgPhosLoader->SetTitle(version);
99   
100   // initialize data members
101   SetDebug(0) ; 
102   fBTE = 0 ; 
103   fPrimaries = 0 ; 
104   fLoadingStatus = "" ; 
105  
106   fESDFileName = rl->GetFileName()  ; // this should be the galice.root file
107   fESDFileName.ReplaceAll("galice.root", "AliESDs.root") ;  
108   fESDFile = 0 ; 
109   fESD = 0 ; 
110   fESDTree = 0 ; 
111   fRawDigits = kFALSE ;
112 }
113
114 //____________________________________________________________________________ 
115 AliPHOSGetter::~AliPHOSGetter()
116 {
117   // dtor
118   if(fgPhosLoader){
119     delete fgPhosLoader ;
120     fgPhosLoader = 0 ;
121   }
122   if(fBTE){
123     delete fBTE ; 
124     fBTE = 0 ;
125   } 
126   if(fPrimaries){
127     fPrimaries->Delete() ; 
128     delete fPrimaries ;
129   } 
130   if (fESD) 
131     delete fESD ; 
132   if (fESDTree) 
133     delete fESDTree ;
134  
135   fgObjGetter = 0;
136 }
137
138 //____________________________________________________________________________ 
139 void AliPHOSGetter::Reset()
140 {
141   // resets things in case the getter is called consecutively with different files
142   // the PHOS Loader is already deleted by the Run Loader
143
144   if (fPrimaries) { 
145     fPrimaries->Delete() ; 
146     delete fPrimaries ;
147   } 
148   fgPhosLoader = 0; 
149   fgObjGetter = 0; 
150 }
151
152 //____________________________________________________________________________ 
153 AliPHOSClusterizer * AliPHOSGetter::Clusterizer()
154
155   // Returns pointer to the Clusterizer task 
156   AliPHOSClusterizer * rv ; 
157   rv =  dynamic_cast<AliPHOSClusterizer *>(PhosLoader()->Reconstructioner()) ;
158   if (!rv) {
159     Event(0, "R") ; 
160     rv =  dynamic_cast<AliPHOSClusterizer*>(PhosLoader()->Reconstructioner()) ;
161   }
162   return rv ; 
163 }
164
165 //____________________________________________________________________________ 
166 TObjArray * AliPHOSGetter::CpvRecPoints() 
167 {
168   // asks the Loader to return the CPV RecPoints container 
169
170   TObjArray * rv = 0 ; 
171   
172   rv = PhosLoader()->CpvRecPoints() ; 
173   if (!rv) {
174     PhosLoader()->MakeRecPointsArray() ;
175     rv = PhosLoader()->CpvRecPoints() ; 
176   }
177   return rv ; 
178 }
179
180 //____________________________________________________________________________ 
181 TClonesArray * AliPHOSGetter::Digits() 
182 {
183   // asks the Loader to return the Digits container 
184
185   TClonesArray * rv = 0 ; 
186   rv = PhosLoader()->Digits() ; 
187
188   if( !rv ) {
189     PhosLoader()->MakeDigitsArray() ; 
190     rv = PhosLoader()->Digits() ;
191   }
192   return rv ; 
193 }
194
195 //____________________________________________________________________________ 
196 AliPHOSDigitizer * AliPHOSGetter::Digitizer() 
197
198   // Returns pointer to the Digitizer task 
199   AliPHOSDigitizer * rv ; 
200   rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
201   if (!rv) {
202     Event(0, "D") ; 
203     rv =  dynamic_cast<AliPHOSDigitizer *>(PhosLoader()->Digitizer()) ;
204   }
205   return rv ; 
206 }
207
208
209 //____________________________________________________________________________ 
210 TObjArray * AliPHOSGetter::EmcRecPoints() 
211 {
212   // asks the Loader to return the EMC RecPoints container 
213
214   TObjArray * rv = 0 ; 
215   
216   rv = PhosLoader()->EmcRecPoints() ; 
217   if (!rv) {
218     PhosLoader()->MakeRecPointsArray() ;
219     rv = PhosLoader()->EmcRecPoints() ; 
220   }
221   return rv ; 
222 }
223
224 //____________________________________________________________________________ 
225 TClonesArray * AliPHOSGetter::TrackSegments() 
226 {
227   // asks the Loader to return the TrackSegments container 
228
229   TClonesArray * rv = 0 ; 
230   
231   rv = PhosLoader()->TrackSegments() ; 
232   if (!rv) {
233     PhosLoader()->MakeTrackSegmentsArray() ;
234     rv = PhosLoader()->TrackSegments() ; 
235   }
236   return rv ; 
237 }
238
239 //____________________________________________________________________________ 
240 AliPHOSTrackSegmentMaker * AliPHOSGetter::TrackSegmentMaker()
241
242   // Returns pointer to the TrackSegmentMaker task 
243   AliPHOSTrackSegmentMaker * rv ; 
244   rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
245   if (!rv) {
246     Event(0, "T") ; 
247     rv =  dynamic_cast<AliPHOSTrackSegmentMaker *>(PhosLoader()->TrackSegmentMaker()) ;
248   }
249   return rv ; 
250 }
251
252 //____________________________________________________________________________ 
253 TClonesArray * AliPHOSGetter::RecParticles() 
254 {
255   // asks the Loader to return the TrackSegments container 
256
257   TClonesArray * rv = 0 ; 
258   
259   rv = PhosLoader()->RecParticles() ; 
260   if (!rv) {
261     PhosLoader()->MakeRecParticlesArray() ;
262     rv = PhosLoader()->RecParticles() ; 
263   }
264   return rv ; 
265 }
266 //____________________________________________________________________________ 
267 void AliPHOSGetter::Event(Int_t event, const char* opt) 
268 {
269   // Reads the content of all Tree's S, D and R
270
271 //   if ( event >= MaxEvent() ) {
272 //     Error("Event", "%d not found in TreeE !", event) ; 
273 //     return ; 
274 //   }
275
276   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
277
278 //   // checks if we are dealing with test-beam data
279 //   TBranch * btb = rl->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
280 //   if(btb){
281 //     if(!fBTE)
282 //       fBTE = new AliPHOSBeamTestEvent() ;
283 //     btb->SetAddress(&fBTE) ;
284 //     btb->GetEntry(event) ;
285 //   }
286 //   else{
287 //     if(fBTE){
288 //       delete fBTE ;
289 //       fBTE = 0 ;
290 //     }
291 //   }
292
293   // Loads the type of object(s) requested
294   
295   rl->GetEvent(event) ;
296
297   if(strstr(opt,"X") || (strcmp(opt,"")==0)){
298     ReadPrimaries() ;
299   }
300   
301   if(strstr(opt,"H")  ){
302     ReadTreeH();
303   }
304   
305   if(strstr(opt,"S")  ){
306     ReadTreeS() ;
307   }
308   
309   if(strstr(opt,"D") ){
310     ReadTreeD() ;
311   }
312   
313   if(strstr(opt,"R") ){
314     ReadTreeR() ;
315   }
316
317   if( strstr(opt,"T") ){
318     ReadTreeT() ;
319   }    
320
321   if( strstr(opt,"P") ){
322     ReadTreeP() ;
323   }    
324
325   if( strstr(opt,"E") ){
326     ReadTreeE(event) ;
327   }
328
329 }
330
331
332 //____________________________________________________________________________ 
333 void AliPHOSGetter::Event(AliRawReader *rawReader, const char* opt) 
334 {
335   // Reads the raw event from rawReader
336
337   if( strstr(opt,"W")  ){
338     ReadRaw(rawReader) ;
339   }    
340  
341 }
342
343
344 //____________________________________________________________________________ 
345 Int_t AliPHOSGetter::EventNumber() const
346   {
347   // return the current event number
348   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
349   return static_cast<Int_t>(rl->GetEventNumber()) ;   
350 }
351
352 //____________________________________________________________________________ 
353   TClonesArray * AliPHOSGetter::Hits()  
354 {
355   // asks the loader to return  the Hits container 
356   
357   TClonesArray * rv = 0 ; 
358   
359   rv = PhosLoader()->Hits() ; 
360   if ( !rv ) {
361     PhosLoader()->LoadHits("read"); 
362     rv = PhosLoader()->Hits() ; 
363   }
364   return rv ; 
365 }
366
367 //____________________________________________________________________________ 
368 AliPHOSGetter * AliPHOSGetter::Instance(const char* alirunFileName, const char* version, Option_t * openingOption) 
369 {
370   // Creates and returns the pointer of the unique instance
371   // Must be called only when the environment has changed
372   
373   if(!fgObjGetter){ // first time the getter is called 
374     fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
375   }
376   else { // the getter has been called previously
377     AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
378     if ( rl->GetFileName() == alirunFileName ) {// the alirunFile has the same name
379       // check if the file is already open
380       TFile * galiceFile = dynamic_cast<TFile *>(gROOT->FindObject(rl->GetFileName()) ) ; 
381       
382       if ( !galiceFile ) 
383         fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;
384       
385       else {  // the file is already open check the version name
386         TString currentVersionName = rl->GetEventFolder()->GetName() ; 
387         TString newVersionName(version) ; 
388         if (currentVersionName == newVersionName) 
389           if(fgDebug)
390             ::Warning( "Instance", "Files with version %s already open", currentVersionName.Data() ) ;  
391         else {
392           fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;      
393         }
394       }
395     }
396     else {
397       AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle()) ; 
398       if ( strstr(version, AliConfig::GetDefaultEventFolderName()) ) // false in case of merging
399         delete rl ; 
400       fgObjGetter = new AliPHOSGetter(alirunFileName, version, openingOption) ;      
401     }
402   }
403   if (!fgObjGetter) 
404     ::Error("AliPHOSGetter::Instance", "Failed to create the PHOS Getter object") ;
405   else 
406     if (fgDebug)
407       Print() ;
408   
409   return fgObjGetter ;
410 }
411
412 //____________________________________________________________________________ 
413 AliPHOSGetter *  AliPHOSGetter::Instance()
414 {
415   // Returns the pointer of the unique instance already defined
416   
417   if(!fgObjGetter && fgDebug)
418      ::Warning("AliPHOSGetter::Instance", "Getter not initialized") ;
419
420    return fgObjGetter ;
421            
422 }
423
424 //____________________________________________________________________________ 
425 Int_t AliPHOSGetter::MaxEvent() const 
426 {
427   // returns the number of events in the run (from TE)
428
429   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
430   return static_cast<Int_t>(rl->GetNumberOfEvents()) ; 
431 }
432
433 //____________________________________________________________________________ 
434 TParticle * AliPHOSGetter::Primary(Int_t index) const
435 {
436   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
437   return rl->Stack()->Particle(index) ; 
438
439
440 //____________________________________________________________________________ 
441 AliPHOS * AliPHOSGetter:: PHOS() const  
442 {
443   // returns the PHOS object 
444   AliPHOSLoader *    loader = 0;
445   static AliPHOSLoader * oldloader = 0;
446   static AliPHOS * phos = 0;
447
448   loader = PhosLoader();
449
450   if ( loader != oldloader) {
451     phos = dynamic_cast<AliPHOS*>(loader->GetModulesFolder()->FindObject("PHOS")) ;
452     oldloader = loader;
453   }
454   if (!phos) 
455     if (fgDebug)
456       Warning("PHOS", "PHOS module not found in module folders: %s", PhosLoader()->GetModulesFolder()->GetName() ) ; 
457   return phos ; 
458 }  
459
460
461
462 //____________________________________________________________________________ 
463 AliPHOSPID * AliPHOSGetter::PID()
464
465   // Returns pointer to the PID task 
466   AliPHOSPID * rv ; 
467   rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
468   if (!rv) {
469     Event(0, "P") ; 
470     rv =  dynamic_cast<AliPHOSPID *>(PhosLoader()->PIDTask()) ;
471   }
472   return rv ; 
473 }
474
475 //____________________________________________________________________________ 
476 AliPHOSGeometry * AliPHOSGetter::PHOSGeometry() const 
477 {
478   // Returns PHOS geometry
479
480   AliPHOSGeometry * rv = 0 ; 
481   if (PHOS() )
482     rv =  PHOS()->GetGeometry() ;
483   return rv ; 
484
485
486 //____________________________________________________________________________ 
487 TClonesArray * AliPHOSGetter::Primaries()  
488 {
489   // creates the Primaries container if needed
490   if ( !fPrimaries ) {
491     if (fgDebug) 
492       Info("Primaries", "Creating a new TClonesArray for primaries") ; 
493     fPrimaries = new TClonesArray("TParticle", 1000) ;
494   } 
495   return fPrimaries ; 
496 }
497
498 //____________________________________________________________________________ 
499 void  AliPHOSGetter::Print() 
500 {
501   // Print usefull information about the getter
502     
503   AliRunLoader * rl = AliRunLoader::GetRunLoader(fgPhosLoader->GetTitle());
504   ::Info( "Print", "gAlice file is %s -- version name is %s", (rl->GetFileName()).Data(), rl->GetEventFolder()->GetName() ) ; 
505 }
506
507 //____________________________________________________________________________ 
508 void AliPHOSGetter::ReadPrimaries()  
509 {
510   // Read Primaries from Kinematics.root
511   
512   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
513   
514   // gets kine tree from the root file (Kinematics.root)
515   if ( ! rl->TreeK() ) { // load treeK the first time
516     rl->LoadKinematics() ;
517   }
518   
519   fNPrimaries = (rl->GetHeader())->GetNtrack(); 
520   if (fgDebug) 
521     Info( "ReadTreeK", "Found %d particles in event # %d", fNPrimaries, EventNumber() ) ; 
522
523
524   // first time creates the container
525   if ( Primaries() ) 
526     fPrimaries->Clear() ; 
527   
528   Int_t index = 0 ; 
529   for (index = 0 ; index < fNPrimaries; index++) { 
530     new ((*fPrimaries)[index]) TParticle(*(Primary(index)));
531   }
532 }
533
534 //____________________________________________________________________________ 
535 Bool_t AliPHOSGetter::OpenESDFile() 
536 {
537   //Open the ESD file    
538   Bool_t rv = kTRUE ; 
539   if (!fESDFile) {
540     fESDFile = TFile::Open(fESDFileName) ;
541     if (!fESDFile ) 
542       return kFALSE ; 
543   }
544   else if (fESDFile->IsOpen()) {
545     fESDFile->Close() ; 
546     fESDFile = TFile::Open(fESDFileName) ;
547   }
548   if (!fESDFile->IsOpen())
549     rv = kFALSE ; 
550   return rv ; 
551 }
552
553 //____________________________________________________________________________ 
554 void AliPHOSGetter::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Int_t & amp, Double_t & time)
555 {
556   // Fits the raw signal time distribution 
557
558   const Int_t kNoiseThreshold = 0 ;
559   Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
560   Double_t signal = 0., signalmax = 0. ;       
561   Double_t energy = time = 0. ; 
562
563   if (lowGainFlag) {
564     timezero1 = timezero2 = signalmax = timemax = 0. ;
565     signalF->FixParameter(0, PHOS()->GetRawFormatLowCharge()) ; 
566     signalF->FixParameter(1, PHOS()->GetRawFormatLowGain()) ; 
567     Int_t index ; 
568     for (index = 0; index < PHOS()->GetRawFormatTimeBins(); index++) {
569       gLowGain->GetPoint(index, time, signal) ; 
570       if (signal > kNoiseThreshold && timezero1 == 0.) 
571         timezero1 = time ;
572       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
573         timezero2 = time ; 
574       if (signal > signalmax) {
575         signalmax = signal ; 
576         timemax   = time ; 
577       }
578     }
579     signalmax /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatLowCharge(), 
580                                                 PHOS()->GetRawFormatLowGain()) ;
581     if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->GetRawFormatTimeMax() * 0.4 ) { // else its noise 
582       signalF->SetParameter(2, signalmax) ; 
583       signalF->SetParameter(3, timezero1) ;         
584       gLowGain->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ; 
585       energy = signalF->GetParameter(2) ; 
586       time   = signalF->GetMaximumX() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
587     }
588   } else {
589     timezero1 = timezero2 = signalmax = timemax = 0. ;
590     signalF->FixParameter(0, PHOS()->GetRawFormatHighCharge()) ; 
591     signalF->FixParameter(1, PHOS()->GetRawFormatHighGain()) ; 
592     Int_t index ; 
593     for (index = 0; index < PHOS()->GetRawFormatTimeBins(); index++) {
594       gHighGain->GetPoint(index, time, signal) ;               
595       if (signal > kNoiseThreshold && timezero1 == 0.) 
596         timezero1 = time ;
597       if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
598         timezero2 = time ; 
599       if (signal > signalmax) {
600         signalmax = signal ;   
601         timemax   = time ; 
602       }
603     }
604     signalmax /= PHOS()->RawResponseFunctionMax(PHOS()->GetRawFormatHighCharge(), 
605                                                 PHOS()->GetRawFormatHighGain()) ;;
606     if ( timezero1 + PHOS()->GetRawFormatTimePeak() < PHOS()->GetRawFormatTimeMax() * 0.4 ) { // else its noise  
607       signalF->SetParameter(2, signalmax) ; 
608       signalF->SetParameter(3, timezero1) ;               
609       gHighGain->Fit(signalF, "QRON", "", 0., timezero2) ; 
610       energy = signalF->GetParameter(2) ; 
611       time   = signalF->GetMaximumX() - PHOS()->GetRawFormatTimePeak() - PHOS()->GetRawFormatTimeTrigger() ;
612     }
613   }
614   
615   if (time == 0. && energy == 0.) 
616     amp = 0 ; 
617   else {
618   AliPHOSDigitizer * digitizer = Digitizer() ; 
619   amp = static_cast<Int_t>( (energy - digitizer->GetEMCpedestal()) / digitizer->GetEMCchannel() + 0.5 ) ; 
620   }
621 }
622
623 //____________________________________________________________________________ 
624 Int_t AliPHOSGetter::ReadRaw(AliRawReader *rawReader)
625 {
626   // reads the raw format data, converts it into digits format and store digits in Digits()
627   // container.
628
629   AliPHOSRawStream in(rawReader);
630  
631   Bool_t first = kTRUE ;
632   
633   TF1 * signalF = new TF1("signal", AliPHOS::RawResponseFunction, 0, PHOS()->GetRawFormatTimeMax(), 4);
634   signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ; 
635
636   Int_t relId[4], id ;
637   Bool_t lowGainFlag = kFALSE ; 
638  
639   TClonesArray * digits = Digits() ;
640   digits->Clear() ; 
641   Int_t idigit = 0 ; 
642   Int_t amp = 0 ; 
643   Double_t time = 0. ; 
644
645   TGraph * gLowGain = new TGraph(PHOS()->GetRawFormatTimeBins()) ; 
646   TGraph * gHighGain= new TGraph(PHOS()->GetRawFormatTimeBins()) ;  
647
648   while ( in.Next() ) { // PHOS entries loop 
649     if ( (in.IsNewRow() || in.IsNewColumn() || in.IsNewModule()) ) {
650       if (!first) {
651         FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
652         if (amp > 0) {
653           new((*digits)[idigit]) AliPHOSDigit( -1, id, amp, time) ;     
654           idigit++ ; 
655         }
656         Int_t index ; 
657         for (index = 0; index < PHOS()->GetRawFormatTimeBins(); index++) {
658           gLowGain->SetPoint(index, index * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 0) ;  
659           gHighGain->SetPoint(index, index * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 0) ; 
660         } 
661       }
662       first = kFALSE ; 
663       relId[0] = in.GetModule() ;
664       if ( relId[0] >= PHOS()->GetRawFormatLowGainOffset() ) { 
665         relId[0] -=  PHOS()->GetRawFormatLowGainOffset() ; 
666         lowGainFlag = kTRUE ; 
667       } else 
668         lowGainFlag = kFALSE ; 
669       relId[1] = 0 ; 
670       relId[2] = in.GetRow() ; 
671       relId[3] = in.GetColumn() ; 
672       PHOSGeometry()->RelToAbsNumbering(relId, id) ;    
673     }
674     if (lowGainFlag)
675       gLowGain->SetPoint(in.GetTime(), 
676                      in.GetTime()* PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 
677                      in.GetSignal()) ;
678     else 
679       gHighGain->SetPoint(in.GetTime(), 
680                          in.GetTime() * PHOS()->GetRawFormatTimeMax() / PHOS()->GetRawFormatTimeBins(), 
681                          in.GetSignal() ) ;
682
683   } // PHOS entries loop
684
685   FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
686   if (amp > 0 ) {
687     new((*digits)[idigit]) AliPHOSDigit( -1, id, amp, time) ;
688     idigit++ ; 
689   }
690   digits->Sort() ; 
691
692   delete signalF ; 
693   delete gLowGain, gHighGain ; 
694   
695   return digits->GetEntriesFast() ; 
696 }
697
698 //   TClonesArray * digits = Digits() ;
699 //   digits->Clear() ; 
700 //   Int_t idigit = 0 ; 
701
702 //   while ( in.Next() ) { // PHOS entries loop 
703  
704 //     Int_t amp = in.GetSignal() ; 
705 //     Double_t time = in.GetTime() ; 
706 //     Int_t relId[4], id ;
707
708 //     relId[0] = in.GetModule() ;
709 //     if ( relId[0] >= PHOS()->GetRawFormatLowGainOffset() ) { 
710 //       relId[0] -=  PHOS()->GetRawFormatLowGainOffset() ; 
711 //     }
712 //     relId[1] = 0 ; 
713 //     relId[2] = in.GetRow() ; 
714 //     relId[3] = in.GetColumn() ; 
715 //     PHOSGeometry()->RelToAbsNumbering(relId, id) ;   
716     
717 //     if (amp > 0 && id >=0 ) {
718 //       new((*digits)[idigit]) AliPHOSDigit( -1, id, amp, time) ;      
719 //       idigit++ ; 
720 //     }
721     
722 //   } // PHOS entries loop
723   
724 //   digits->Sort() ; 
725   
726 //   return digits->GetEntriesFast() ; 
727 // }
728 //____________________________________________________________________________ 
729 Int_t AliPHOSGetter::ReadTreeD()
730 {
731   // Read the Digits
732   
733   PhosLoader()->CleanDigits() ;    
734   PhosLoader()->LoadDigits("UPDATE") ;
735   PhosLoader()->LoadDigitizer("UPDATE") ;
736   return Digits()->GetEntries() ; 
737 }
738
739 //____________________________________________________________________________ 
740 Int_t AliPHOSGetter::ReadTreeH()
741 {
742   // Read the Hits
743   PhosLoader()->CleanHits() ;
744   // gets TreeH from the root file (PHOS.Hit.root)
745   //if ( !IsLoaded("H") ) {
746     PhosLoader()->LoadHits("UPDATE") ;
747   //  SetLoaded("H") ; 
748   //}  
749   return Hits()->GetEntries() ; 
750 }
751
752 //____________________________________________________________________________ 
753 Int_t AliPHOSGetter::ReadTreeR()
754 {
755   // Read the RecPoints
756   
757   PhosLoader()->CleanRecPoints() ;
758   // gets TreeR from the root file (PHOS.RecPoints.root)
759   //if ( !IsLoaded("R") ) {
760     PhosLoader()->LoadRecPoints("UPDATE") ;
761     PhosLoader()->LoadClusterizer("UPDATE") ;
762     //  SetLoaded("R") ; 
763     //}
764
765   return EmcRecPoints()->GetEntries() ; 
766 }
767
768 //____________________________________________________________________________ 
769 Int_t AliPHOSGetter::ReadTreeT()
770 {
771   // Read the TrackSegments
772   
773   PhosLoader()->CleanTracks() ; 
774   // gets TreeT from the root file (PHOS.TrackSegments.root)
775   //if ( !IsLoaded("T") ) {
776     PhosLoader()->LoadTracks("UPDATE") ;
777     PhosLoader()->LoadTrackSegmentMaker("UPDATE") ;
778     //    SetLoaded("T") ; 
779     //}
780
781   return TrackSegments()->GetEntries() ; 
782 }
783 //____________________________________________________________________________ 
784 Int_t AliPHOSGetter::ReadTreeP()
785 {
786   // Read the RecParticles
787   
788   PhosLoader()->CleanRecParticles() ; 
789
790   // gets TreeT from the root file (PHOS.TrackSegments.root)
791   //  if ( !IsLoaded("P") ) {
792     PhosLoader()->LoadRecParticles("UPDATE") ;
793     PhosLoader()->LoadPID("UPDATE") ;
794     //  SetLoaded("P") ; 
795     //}
796
797   return RecParticles()->GetEntries() ; 
798 }
799 //____________________________________________________________________________ 
800 Int_t AliPHOSGetter::ReadTreeS()
801 {
802   // Read the SDigits
803   
804   PhosLoader()->CleanSDigits() ; 
805   // gets TreeS from the root file (PHOS.SDigits.root)
806   //if ( !IsLoaded("S") ) {
807     PhosLoader()->LoadSDigits("READ") ;
808     PhosLoader()->LoadSDigitizer("READ") ;
809     //  SetLoaded("S") ; 
810     //}
811
812   return SDigits()->GetEntries() ; 
813 }
814
815 //____________________________________________________________________________ 
816 Int_t AliPHOSGetter::ReadTreeE(Int_t event)
817 {
818   // Read the ESD
819   
820   // gets esdTree from the root file (AliESDs.root)
821   if (!fESDFile)
822     if ( !OpenESDFile() ) 
823       return -1 ; 
824
825   fESDTree = static_cast<TTree*>(fESDFile->Get("esdTree")) ; 
826   fESD = new AliESD;
827    if (!fESDTree) {
828
829      Error("ReadTreeE", "no ESD tree found");
830      return -1;
831    }
832    fESDTree->SetBranchAddress("ESD", &fESD);
833    fESDTree->GetEvent(event);
834
835    return event ; 
836 }
837
838 //____________________________________________________________________________ 
839 TClonesArray * AliPHOSGetter::SDigits() 
840 {
841   // asks the Loader to return the Digits container 
842
843   TClonesArray * rv = 0 ; 
844   
845   rv = PhosLoader()->SDigits() ; 
846   if (!rv) {
847     PhosLoader()->MakeSDigitsArray() ;
848     rv = PhosLoader()->SDigits() ; 
849   }
850   return rv ; 
851 }
852
853 //____________________________________________________________________________ 
854 AliPHOSSDigitizer * AliPHOSGetter::SDigitizer()
855
856   // Returns pointer to the SDigitizer task 
857   AliPHOSSDigitizer * rv ; 
858   rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
859   if (!rv) {
860     Event(0, "S") ; 
861     rv =  dynamic_cast<AliPHOSSDigitizer *>(PhosLoader()->SDigitizer()) ;
862   }
863   return rv ; 
864 }
865
866 //____________________________________________________________________________ 
867 TParticle * AliPHOSGetter::Secondary(const TParticle* p, Int_t index) const
868 {
869   // Return first (index=1) or second (index=2) secondary particle of primary particle p 
870
871   if(index <= 0) 
872     return 0 ;
873   if(index > 2)
874     return 0 ;
875
876   if(p) {
877   Int_t daughterIndex = p->GetDaughter(index-1) ; 
878   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
879   return  rl->GetAliRun()->GetMCApp()->Particle(daughterIndex) ; 
880   }
881   else
882     return 0 ;
883 }
884
885 //____________________________________________________________________________ 
886 void AliPHOSGetter::Track(Int_t itrack) 
887 {
888   // Read the first entry of PHOS branch in hit tree gAlice->TreeH()
889  
890  AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
891
892   if( !TreeH() ) // load treeH the first time
893     rl->LoadHits() ;
894
895   // first time create the container
896   TClonesArray * hits = Hits() ; 
897   if ( hits ) 
898     hits->Clear() ; 
899
900   TBranch * phosbranch = dynamic_cast<TBranch*>(TreeH()->GetBranch("PHOS")) ; 
901   phosbranch->SetAddress(&hits) ;
902   phosbranch->GetEntry(itrack) ;
903 }
904
905 //____________________________________________________________________________ 
906 TTree * AliPHOSGetter::TreeD() const 
907 {
908   // Returns pointer to the Digits Tree
909   TTree * rv = 0 ; 
910   rv = PhosLoader()->TreeD() ; 
911   if ( !rv ) {
912     PhosLoader()->MakeTree("D");
913     rv = PhosLoader()->TreeD() ;
914   } 
915   
916   return rv ; 
917 }
918
919 //____________________________________________________________________________ 
920 TTree * AliPHOSGetter::TreeH() const 
921 {
922   // Returns pointer to the Hits Tree
923   TTree * rv = 0 ; 
924   rv = PhosLoader()->TreeH() ; 
925   if ( !rv ) {
926     PhosLoader()->MakeTree("H");
927     rv = PhosLoader()->TreeH() ;
928   } 
929   
930   return rv ; 
931 }
932
933 //____________________________________________________________________________ 
934 TTree * AliPHOSGetter::TreeR() const 
935 {
936   // Returns pointer to the RecPoints Tree
937   TTree * rv = 0 ; 
938   rv = PhosLoader()->TreeR() ; 
939   if ( !rv ) {
940     PhosLoader()->MakeTree("R");
941     rv = PhosLoader()->TreeR() ;
942   } 
943   
944   return rv ; 
945 }
946
947 //____________________________________________________________________________ 
948 TTree * AliPHOSGetter::TreeT() const 
949 {
950   // Returns pointer to the TrackSegments Tree
951   TTree * rv = 0 ; 
952   rv = PhosLoader()->TreeT() ; 
953   if ( !rv ) {
954     PhosLoader()->MakeTree("T");
955     rv = PhosLoader()->TreeT() ;
956   } 
957   
958   return rv ; 
959 }
960 //____________________________________________________________________________ 
961 TTree * AliPHOSGetter::TreeP() const 
962 {
963   // Returns pointer to the RecParticles  Tree
964   TTree * rv = 0 ; 
965   rv = PhosLoader()->TreeP() ; 
966   if ( !rv ) {
967     PhosLoader()->MakeTree("P");
968     rv = PhosLoader()->TreeP() ;
969   } 
970   
971   return rv ; 
972 }
973
974 //____________________________________________________________________________ 
975 TTree * AliPHOSGetter::TreeS() const 
976
977  // Returns pointer to the SDigits Tree
978   TTree * rv = 0 ; 
979   rv = PhosLoader()->TreeS() ; 
980   if ( !rv ) {
981     PhosLoader()->MakeTree("S");
982     rv = PhosLoader()->TreeS() ;
983   } 
984   
985   return rv ; 
986 }
987
988 //____________________________________________________________________________ 
989 Bool_t AliPHOSGetter::VersionExists(TString & opt) const
990 {
991   // checks if the version with the present name already exists in the same directory
992
993   Bool_t rv = kFALSE ;
994  
995   AliRunLoader * rl = AliRunLoader::GetRunLoader(PhosLoader()->GetTitle());
996   TString version( rl->GetEventFolder()->GetName() ) ; 
997
998   opt.ToLower() ; 
999   
1000   if ( opt == "sdigits") {
1001     // add the version name to the root file name
1002     TString fileName( PhosLoader()->GetSDigitsFileName() ) ; 
1003     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
1004       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1005     if ( !(gSystem->AccessPathName(fileName)) ) { 
1006       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
1007       rv = kTRUE ; 
1008     }
1009     PhosLoader()->SetSDigitsFileName(fileName) ;
1010   }
1011
1012   if ( opt == "digits") {
1013     // add the version name to the root file name
1014     TString fileName( PhosLoader()->GetDigitsFileName() ) ; 
1015     if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
1016       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
1017     if ( !(gSystem->AccessPathName(fileName)) ) {
1018       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;  
1019       rv = kTRUE ; 
1020     }
1021   }
1022
1023   return rv ;
1024
1025 }
1026
1027 //____________________________________________________________________________ 
1028 UShort_t AliPHOSGetter::EventPattern(void) const
1029 {
1030   // Return the pattern (trigger bit register) of the beam-test event
1031   if(fBTE)
1032     return fBTE->GetPattern() ;
1033   else
1034     return 0 ;
1035 }
1036 //____________________________________________________________________________ 
1037 Float_t AliPHOSGetter::BeamEnergy(void) const
1038 {
1039   // Return the beam energy of the beam-test event
1040   if(fBTE)
1041     return fBTE->GetBeamEnergy() ;
1042   else
1043     return 0 ;
1044 }