restore threshold getters after parameter dynamics update (fw v. >= A012)
[u/mrichter/AliRoot.git] / FMD / AliFMDInput.cxx
1 /**************************************************************************
2  * Copyright(c) 2004, 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 /* $Id$ */
16 /** @file    AliFMDInput.cxx
17     @author  Christian Holm Christensen <cholm@nbi.dk>
18     @date    Mon Mar 27 12:42:40 2006
19     @brief   FMD utility classes for reading FMD data
20 */
21 //___________________________________________________________________
22 //
23 // The classes defined here, are utility classes for reading in data
24 // for the FMD.  They are  put in a seperate library to not polute the
25 // normal libraries.  The classes are intended to be used as base
26 // classes for customized class that do some sort of analysis on the
27 // various types of data produced by the FMD. 
28 //
29 // Latest changes by Christian Holm Christensen
30 //
31 #include "AliFMDInput.h"        // ALIFMDHIT_H
32 #include "AliFMDDebug.h"        // ALIFMDDEBUG_H ALILOG_H
33 #include "AliLoader.h"          // ALILOADER_H
34 #include "AliRunLoader.h"       // ALIRUNLOADER_H
35 #include "AliRun.h"             // ALIRUN_H
36 #include "AliStack.h"           // ALISTACK_H
37 #include "AliRawReaderFile.h"   // ALIRAWREADERFILE_H
38 #include "AliRawReaderRoot.h"   // ALIRAWREADERROOT_H
39 #include "AliRawReaderDate.h"   // ALIRAWREADERDATE_H
40 #include "AliRawEventHeaderBase.h" 
41 #include "AliFMD.h"             // ALIFMD_H
42 #include "AliFMDHit.h"          // ALIFMDHIT_H
43 #include "AliFMDDigit.h"        // ALIFMDDigit_H
44 #include "AliFMDSDigit.h"       // ALIFMDDigit_H
45 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
46 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
47 #include "AliFMDGeometry.h"
48 #include <AliESD.h>
49 #include <AliESDFMD.h>
50 #include <AliESDEvent.h>
51 #include <AliCDBManager.h>
52 #include <AliCDBEntry.h>
53 #include <AliAlignObjParams.h>
54 #include <AliTrackReference.h>
55 #include <TTree.h>              // ROOT_TTree
56 #include <TChain.h>             // ROOT_TChain
57 #include <TParticle.h>          // ROOT_TParticle
58 #include <TString.h>            // ROOT_TString
59 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
60 #include <TMath.h>              // ROOT_TMath
61 #include <TGeoManager.h>        // ROOT_TGeoManager 
62 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
63 #include <Riostream.h>          // ROOT_Riostream
64 #include <TFile.h>              // ROOT_TFile
65 #include <TStreamerInfo.h>
66 #include <TArrayF.h>
67
68 //____________________________________________________________________
69 ClassImp(AliFMDInput)
70 #if 0
71   ; // This is here to keep Emacs for indenting the next line
72 #endif
73
74 //____________________________________________________________________
75 const AliFMDInput::ETrees AliFMDInput::fgkAllLoads[] = { kHits,       
76                                                          kKinematics,
77                                                          kDigits,    
78                                                          kSDigits,   
79                                                          kHeader,    
80                                                          kRecPoints, 
81                                                          kESD,       
82                                                          kRaw,       
83                                                          kGeometry,  
84                                                          kTrackRefs, 
85                                                          kRawCalib,  
86                                                          kUser };
87
88 //____________________________________________________________________
89 AliFMDInput::AliFMDInput()
90   : TNamed("AliFMDInput", "Input handler for various FMD data"), 
91     fGAliceFile(""), 
92     fLoader(0),
93     fRun(0), 
94     fStack(0),
95     fFMDLoader(0), 
96     fReader(0),
97     fFMDReader(0),
98     fFMD(0),
99     fESD(0),
100     fESDEvent(0),
101     fTreeE(0),
102     fTreeH(0),
103     fTreeTR(0),
104     fTreeD(0),
105     fTreeS(0),
106     fTreeR(0), 
107     fTreeA(0), 
108     fChainE(0),
109     fArrayE(0),
110     fArrayH(0),
111     fArrayTR(0), 
112     fArrayD(0),
113     fArrayS(0), 
114     fArrayR(0), 
115     fArrayA(0), 
116     fHeader(0),
117     fGeoManager(0),
118     fTreeMask(0), 
119     fRawFile(""),
120     fInputDir("."),
121     fIsInit(kFALSE),
122     fEventCount(0), 
123     fNEvents(-1)
124 {
125
126   // Constructor of an FMD input object.  Specify what data to read in
127   // using the AddLoad member function.   Sub-classes should at a
128   // minimum overload the member function Event.   A full job can be
129   // executed using the member function Run. 
130 }
131
132   
133
134 //____________________________________________________________________
135 AliFMDInput::AliFMDInput(const char* gAliceFile)
136   : TNamed("AliFMDInput", "Input handler for various FMD data"), 
137     fGAliceFile(gAliceFile),
138     fLoader(0),
139     fRun(0), 
140     fStack(0),
141     fFMDLoader(0), 
142     fReader(0),
143     fFMDReader(0),
144     fFMD(0),
145     fESD(0),
146     fESDEvent(0),
147     fTreeE(0),
148     fTreeH(0),
149     fTreeTR(0),
150     fTreeD(0),
151     fTreeS(0),
152     fTreeR(0), 
153     fTreeA(0), 
154     fChainE(0),
155     fArrayE(0),
156     fArrayH(0),
157     fArrayTR(0), 
158     fArrayD(0),
159     fArrayS(0), 
160     fArrayR(0), 
161     fArrayA(0), 
162     fHeader(0),
163     fGeoManager(0),
164     fTreeMask(0), 
165     fRawFile(""),
166     fInputDir("."),
167     fIsInit(kFALSE),
168     fEventCount(0),
169     fNEvents(-1)
170 {
171   
172   // Constructor of an FMD input object.  Specify what data to read in
173   // using the AddLoad member function.   Sub-classes should at a
174   // minimum overload the member function Event.   A full job can be
175   // executed using the member function Run. 
176 }
177
178 //____________________________________________________________________
179 void
180 AliFMDInput::SetLoads(UInt_t mask)
181 {
182   for (UInt_t i = 0; i < sizeof(mask); i++) { 
183     if (!(mask & (1 << i))) continue;
184     const ETrees *ptype = fgkAllLoads;
185     do { 
186       ETrees type = *ptype;
187       if (i != UInt_t(type)) continue;
188       AddLoad(type);
189       break;
190     } while (*ptype++ != kUser);
191   }
192 }
193      
194 //____________________________________________________________________
195 void
196 AliFMDInput::SetLoads(const char* what) 
197 {
198   TString    l(what);
199   TObjArray* ll = l.Tokenize(", ");
200   TIter      next(ll);
201   TObject*   os = 0;
202   while ((os = next())) { 
203     ETrees type = ParseLoad(os->GetName());
204     AddLoad(type);
205   }
206 }
207     
208
209 //____________________________________________________________________
210 AliFMDInput::ETrees
211 AliFMDInput::ParseLoad(const char* what)
212 {
213   TString opt(what);
214   opt.ToLower();
215   const ETrees* ptype = fgkAllLoads;
216   do { 
217     ETrees  type = *ptype;
218     if (opt.Contains(TreeName(type,true), TString::kIgnoreCase)) 
219       return type;
220   } while (*ptype++ != kUser);
221   return kUser;
222 }
223 //____________________________________________________________________
224 const char*
225 AliFMDInput::LoadedString(Bool_t dataOnly) const
226 {
227   static TString ret;
228   if (!ret.IsNull()) return ret.Data();
229
230   const ETrees* ptype = fgkAllLoads;
231   do { 
232     ETrees type = *ptype;
233     if (dataOnly && 
234         (type == kKinematics || 
235          type == kHeader     || 
236          type == kGeometry   || 
237          type == kTrackRefs)) continue;
238     if (!IsLoaded(*ptype)) continue;
239     
240     if (!ret.IsNull()) ret.Append(",");
241     ret.Append(TreeName(type));
242   } while (*ptype++ != kUser);
243   return ret.Data();
244 }
245
246 //____________________________________________________________________
247 const char* 
248 AliFMDInput::TreeName(ETrees tree, Bool_t shortest) 
249 {
250   if (shortest) {
251     switch (tree) {
252     case kHits:         return "hit";        
253     case kKinematics:   return "kin";
254     case kDigits:       return "dig";     
255     case kSDigits:      return "sdig";    
256     case kHeader:       return "hea";     
257     case kRecPoints:    return "recp";  
258     case kESD:          return "esd";        
259     case kRaw:          return "raw";        
260     case kGeometry:     return "geo";   
261     case kTrackRefs:    return "trackr";  
262     case kRawCalib:     return "rawc";   
263     case kUser:         return "user"; 
264     }
265     return 0;
266   }
267   switch (tree) {
268   case kHits:           return "Hits";        
269   case kKinematics:     return "Kinematics"; 
270   case kDigits:         return "Digits";     
271   case kSDigits:        return "SDigits";    
272   case kHeader:         return "Header";     
273   case kRecPoints:      return "RecPoints";  
274   case kESD:            return "ESD";        
275   case kRaw:            return "Raw";        
276   case kGeometry:       return "Geometry";   
277   case kTrackRefs:      return "TrackRefs";  
278   case kRawCalib:       return "RawCalib";   
279   case kUser:           return "User"; 
280   }
281   return 0;
282 }
283
284 //____________________________________________________________________
285 Int_t
286 AliFMDInput::NEvents() const 
287 {
288   // Get number of events
289   if (IsLoaded(kRaw) || 
290       IsLoaded(kRawCalib)) return fReader->GetNumberOfEvents();
291   if (fChainE) return fChainE->GetEntriesFast();
292   if (fTreeE) return fTreeE->GetEntries();
293   return -1;
294 }
295
296 //____________________________________________________________________
297 Bool_t
298 AliFMDInput::Init()
299 {
300   // Initialize the object.  Get the needed loaders, and such. 
301
302   // Check if we have been initialized
303   if (fIsInit) { 
304     AliWarning("Already initialized");
305     return fIsInit;
306   }
307   TString what;
308   const ETrees* ptype = fgkAllLoads;
309   do { 
310     ETrees type = *ptype;
311     what.Append(Form("\n\t%-20s: %s", TreeName(type), 
312                      IsLoaded(type) ? "yes" : "no"));
313   } while (*ptype++ != kUser);
314   
315   Info("Init","Initialising w/mask 0x%04x%s", fTreeMask, what.Data());
316   // Get the run 
317   if (IsLoaded(kDigits)     ||
318       IsLoaded(kSDigits)    || 
319       IsLoaded(kKinematics) || 
320       IsLoaded(kTrackRefs)  || 
321       IsLoaded(kHeader)) {
322     if (!gSystem->FindFile(".:/", fGAliceFile)) {
323       AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
324     }
325     else {
326       fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
327       if (!fLoader) {
328         AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
329         return kFALSE;
330       }
331       AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
332
333       if  (fLoader->LoadgAlice()) return kFALSE;
334       
335       fRun = fLoader->GetAliRun();
336       
337       // Get the FMD 
338       fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
339       if (!fFMD) {
340         AliError("Failed to get detector FMD from loader");
341         return kFALSE;
342       }
343       
344       // Get the FMD loader
345       fFMDLoader = fLoader->GetLoader("FMDLoader");
346       if (!fFMDLoader) {
347         AliError("Failed to get detector FMD loader from loader");
348         return kFALSE;
349       }
350       if (fLoader->LoadHeader()) { 
351         AliError("Failed to get event header information from loader");
352         return kFALSE;
353       }
354       fTreeE = fLoader->TreeE();
355     }
356   }
357
358   // Optionally, get the ESD files
359   if (IsLoaded(kESD)) {
360     fChainE = MakeChain("ESD", fInputDir, true);
361     fESDEvent = new AliESDEvent();
362     fESDEvent->ReadFromTree(fChainE);
363     //    fChainE->SetBranchAddress("ESD", &fMainESD);
364     
365   }
366     
367   if (IsLoaded(kRaw) || 
368       IsLoaded(kRawCalib)) {
369     AliInfo("Getting FMD raw data digits");
370     fArrayA = new TClonesArray("AliFMDDigit");
371 #if 0
372     if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
373       fReader = new AliRawReaderRoot(fRawFile.Data());
374     else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
375       fReader = new AliRawReaderDate(fRawFile.Data());
376     else
377       fReader = new AliRawReaderFile(-1);
378 #else
379     if(!fRawFile.IsNull()) 
380       fReader = AliRawReader::Create(fRawFile.Data());
381     else 
382       fReader = new AliRawReaderFile(-1);
383 #endif
384     fFMDReader = new AliFMDRawReader(fReader, 0);
385   }
386   
387   // Optionally, get the geometry 
388   if (IsLoaded(kGeometry)) {
389     TString fname;
390     if (fRun) {
391       fname = gSystem->DirName(fGAliceFile);
392       fname.Append("/geometry.root");
393     }
394     if (!gSystem->AccessPathName(fname.Data())) 
395       fname = "";
396     AliCDBManager* cdb   = AliCDBManager::Instance();
397     if (!cdb->IsDefaultStorageSet()) {
398       cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
399       cdb->SetRun(0);
400     }
401
402     AliGeomManager::LoadGeometry(fname.IsNull() ? 0 : fname.Data());
403
404     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
405     if (align) {
406       AliInfo("Got alignment data from CDB");
407       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
408       if (!array) {
409         AliWarning("Invalid align data from CDB");
410       }
411       else {
412         Int_t nAlign = array->GetEntries();
413         for (Int_t i = 0; i < nAlign; i++) {
414           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
415           if (!a->ApplyToGeometry()) {
416             AliWarning(Form("Failed to apply alignment to %s", 
417                             a->GetSymName()));
418           }
419         }
420       }
421     }
422     AliFMDGeometry* geom = AliFMDGeometry::Instance();
423     geom->Init();
424     geom->InitTransformations();
425   }
426
427   fEventCount = 0;
428   fIsInit = kTRUE;
429   return fIsInit;
430 }
431
432 //____________________________________________________________________
433 Bool_t
434 AliFMDInput::Begin(Int_t event)
435 {
436   // Called at the begining of each event.  Per default, it gets the
437   // data trees and gets pointers to the output arrays.   Users can
438   // overload this, but should call this member function in the
439   // overloaded member function of the derived class. 
440
441   // Check if we have been initialized
442   if (!fIsInit) { 
443     AliError("Not initialized");
444     return fIsInit;
445   }
446
447   // Get the event 
448   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
449
450   // Possibly load global kinematics information 
451   if (IsLoaded(kKinematics)) {
452     // AliInfo("Getting kinematics");
453     if (fLoader->LoadKinematics("READ")) return kFALSE;
454     fStack = fLoader->Stack();
455   }
456
457   // Possibly load FMD Hit information 
458   if (IsLoaded(kHits)) {
459     // AliInfo("Getting FMD hits");
460     if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
461     fTreeH = fFMDLoader->TreeH();
462     if (!fArrayH) fArrayH = fFMD->Hits(); 
463   }
464   
465   // Possibly load FMD TrackReference information 
466   if (IsLoaded(kTrackRefs)) {
467     // AliInfo("Getting FMD hits");
468     if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
469     fTreeTR = fLoader->TreeTR();
470     if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
471     fTreeTR->SetBranchAddress("TrackReferences",  &fArrayTR);
472   }
473   
474   // Possibly load heaedr information 
475   if (IsLoaded(kHeader)) {
476     // AliInfo("Getting FMD hits");
477     if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
478     fHeader = fLoader->GetHeader();
479   }
480
481   // Possibly load FMD Digit information 
482   if (IsLoaded(kDigits)) {
483     // AliInfo("Getting FMD digits");
484     if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
485     fTreeD = fFMDLoader->TreeD();
486     if (fTreeD) {
487       if (!fArrayD) fArrayD = fFMD->Digits();
488     }
489     else {
490       fArrayD = 0;
491       AliWarning(Form("Failed to load FMD Digits"));
492     } 
493   }
494
495   // Possibly load FMD Sdigit information 
496   if (IsLoaded(kSDigits)) {
497     // AliInfo("Getting FMD summable digits");
498     if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) { 
499       AliWarning("Failed to load SDigits!");
500       return kFALSE;
501     }
502     fTreeS = fFMDLoader->TreeS();
503     if (!fArrayS) fArrayS = fFMD->SDigits();
504   }
505
506   // Possibly load FMD RecPoints information 
507   if (IsLoaded(kRecPoints)) {
508     // AliInfo("Getting FMD reconstructed points");
509     if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
510     fTreeR = fFMDLoader->TreeR();
511     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
512     fTreeR->SetBranchAddress("FMD",  &fArrayR);
513   }  
514
515   // Possibly load FMD ESD information 
516   if (IsLoaded(kESD)) {
517     // AliInfo("Getting FMD event summary data");
518     Int_t read = fChainE->GetEntry(event);
519     if (read <= 0) return kFALSE;
520     fESD = fESDEvent->GetFMDData();
521     if (!fESD) return kFALSE;
522   }
523
524   // Possibly load FMD Digit information 
525   if (IsLoaded(kRaw) || IsLoaded(kRawCalib)) {
526     Bool_t mon = fRawFile.Contains("mem://");
527     // AliInfo("Getting FMD raw data digits");
528     if (mon) std::cout << "Waiting for event ..." << std::flush;
529     do { 
530       if (!fReader->NextEvent()) { 
531         if (mon) { 
532           gSystem->Sleep(3);
533           continue;
534         }
535         return kFALSE;
536       }
537       UInt_t eventType = fReader->GetType();
538       if(eventType == AliRawEventHeaderBase::kPhysicsEvent ||
539          eventType == AliRawEventHeaderBase::kCalibrationEvent) 
540         break;
541     } while (true);
542     if (mon) std::cout << "got it" << std::endl;
543     // AliFMDRawReader r(fReader, 0);
544     fArrayA->Clear();
545     fFMDReader->ReadAdcs(fArrayA);
546     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
547   }
548   fEventCount++;
549   return kTRUE;
550 }
551
552
553 //____________________________________________________________________
554 Bool_t 
555 AliFMDInput::Event()
556 {
557   // Process one event.  The default implementation one or more of 
558   //
559   //   -  ProcessHits       if the hits are loaded. 
560   //   -  ProcessDigits     if the digits are loaded. 
561   //   -  ProcessSDigits    if the sumbable digits are loaded. 
562   //   -  ProcessRecPoints  if the reconstructed points are loaded. 
563   //   -  ProcessESD        if the event summary data is loaded
564   // 
565   if (IsLoaded(kHits))     if (!ProcessHits())      return kFALSE; 
566   if (IsLoaded(kTrackRefs))if (!ProcessTrackRefs()) return kFALSE; 
567   if (IsLoaded(kKinematics) && 
568       IsLoaded(kHits))     if (!ProcessTracks())    return kFALSE; 
569   if (IsLoaded(kKinematics))if (!ProcessStack())     return kFALSE; 
570   if (IsLoaded(kSDigits))  if (!ProcessSDigits())   return kFALSE;
571   if (IsLoaded(kDigits))   if (!ProcessDigits())    return kFALSE;
572   if (IsLoaded(kRaw))      if (!ProcessRawDigits()) return kFALSE;
573   if (IsLoaded(kRawCalib)) if (!ProcessRawCalibDigits())return kFALSE;
574   if (IsLoaded(kRecPoints))if (!ProcessRecPoints()) return kFALSE;
575   if (IsLoaded(kESD))      if (!ProcessESDs())      return kFALSE;
576   if (IsLoaded(kUser))     if (!ProcessUsers())     return kFALSE;
577   
578   return kTRUE;
579 }
580
581 //____________________________________________________________________
582 Bool_t 
583 AliFMDInput::ProcessHits()
584 {
585   // Read the hit tree, and pass each hit to the member function
586   // ProcessHit.
587   if (!fTreeH) {
588     AliError("No hit tree defined");
589     return kFALSE;
590   }
591   if (!fArrayH) {
592     AliError("No hit array defined");
593     return kFALSE;
594   }
595
596   Int_t nTracks = fTreeH->GetEntries();
597   for (Int_t i = 0; i < nTracks; i++) {
598     Int_t hitRead  = fTreeH->GetEntry(i);
599     if (hitRead <= 0) continue;
600
601     Int_t nHit = fArrayH->GetEntries();
602     if (nHit <= 0) continue;
603   
604     for (Int_t j = 0; j < nHit; j++) {
605       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
606       if (!hit) continue;
607
608       TParticle* track = 0;
609       if (IsLoaded(kKinematics) && fStack) {
610         Int_t trackno = hit->Track();
611         track = fStack->Particle(trackno);
612       }
613       if (!ProcessHit(hit, track)) return kFALSE;
614     }    
615   }
616   return kTRUE;
617 }
618 //____________________________________________________________________
619 Bool_t 
620 AliFMDInput::ProcessTrackRefs()
621 {
622   // Read the reconstrcted points tree, and pass each reconstruction
623   // object (AliFMDRecPoint) to either ProcessRecPoint.
624   if (!fTreeTR) {
625     AliError("No track reference tree defined");
626     return kFALSE;
627   }
628   if (!fArrayTR) {
629     AliError("No track reference array defined");
630     return kFALSE;
631   }
632
633   Int_t nEv = fTreeTR->GetEntries();
634   for (Int_t i = 0; i < nEv; i++) {
635     Int_t trRead  = fTreeTR->GetEntry(i);
636     if (trRead <= 0) continue;
637     Int_t nTrackRefs = fArrayTR->GetEntries();
638     for (Int_t j = 0; j < nTrackRefs; j++) {
639       AliTrackReference* trackRef = 
640         static_cast<AliTrackReference*>(fArrayTR->At(j));
641       if (!trackRef) continue;
642       // if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
643       TParticle* track = 0;
644       if (IsLoaded(kKinematics) && fStack) {
645         Int_t trackno = trackRef->GetTrack();
646         track = fStack->Particle(trackno);
647       }
648       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
649     }    
650   }
651   return kTRUE;
652 }
653 //____________________________________________________________________
654 Bool_t 
655 AliFMDInput::ProcessTracks()
656 {
657   // Read the hit tree, and pass each hit to the member function
658   // ProcessTrack.
659   if (!fStack) {
660     AliError("No track tree defined");
661     return kFALSE;
662   }
663   if (!fTreeH) {
664     AliError("No hit tree defined");
665     return kFALSE;
666   }
667   if (!fArrayH) {
668     AliError("No hit array defined");
669     return kFALSE;
670   }
671
672   // Int_t nTracks = fStack->GetNtrack();
673   Int_t nTracks = fTreeH->GetEntries();
674   for (Int_t i = 0; i < nTracks; i++) {
675     Int_t      trackno = nTracks - i - 1;
676     TParticle* track   = fStack->Particle(trackno);
677     if (!track) continue;
678
679     // Get the hits for this track. 
680     Int_t hitRead  = fTreeH->GetEntry(i);
681     Int_t nHit     = fArrayH->GetEntries();
682     if (nHit == 0 || hitRead <= 0) { 
683       // Let user code see the track, even if there's no hits. 
684       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
685       continue;
686     }
687
688     // Loop over the hits corresponding to this track.
689     for (Int_t j = 0; j < nHit; j++) {
690       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
691       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
692     }   
693   }
694   return kTRUE;
695 }
696 //____________________________________________________________________
697 Bool_t 
698 AliFMDInput::ProcessStack()
699 {
700   // Read the hit tree, and pass each hit to the member function
701   // ProcessTrack.
702   if (!fStack) {
703     AliError("No track tree defined");
704     return kFALSE;
705   }
706   Int_t nTracks = fStack->GetNtrack();
707   for (Int_t i = 0; i < nTracks; i++) {
708     Int_t      trackno = nTracks - i - 1;
709     TParticle* track   = fStack->Particle(trackno);
710     if (!track) continue;
711
712     if (!ProcessParticle(trackno, track)) return kFALSE;
713   }
714   return kTRUE;
715 }
716 //____________________________________________________________________
717 Bool_t 
718 AliFMDInput::ProcessDigits()
719 {
720   // Read the digit tree, and pass each digit to the member function
721   // ProcessDigit.
722   if (!fTreeD) {
723     AliError("No digit tree defined");
724     return kFALSE;
725   }
726   if (!fArrayD) {
727     AliError("No digit array defined");
728     return kFALSE;
729   }
730
731   Int_t nEv = fTreeD->GetEntries();
732   for (Int_t i = 0; i < nEv; i++) {
733     Int_t digitRead  = fTreeD->GetEntry(i);
734     if (digitRead <= 0) continue;
735     Int_t nDigit = fArrayD->GetEntries();
736     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
737     if (nDigit <= 0) continue;
738     for (Int_t j = 0; j < nDigit; j++) {
739       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
740       if (!digit) continue;
741       if (!ProcessDigit(digit)) return kFALSE;
742     }    
743   }
744   return kTRUE;
745 }
746
747 //____________________________________________________________________
748 Bool_t 
749 AliFMDInput::ProcessSDigits()
750 {
751   // Read the summable digit tree, and pass each sumable digit to the
752   // member function ProcessSdigit.
753   if (!fTreeS) {
754     AliWarning("No sdigit tree defined");
755     return kTRUE; // Empty SDigits is fine
756   }
757   if (!fArrayS) {
758     AliWarning("No sdigit array defined");
759     return kTRUE; // Empty SDigits is fine
760   }
761
762   Int_t nEv = fTreeS->GetEntries();
763   for (Int_t i = 0; i < nEv; i++) {
764     Int_t sdigitRead  = fTreeS->GetEntry(i);
765     if (sdigitRead <= 0) { 
766       AliInfo(Form("Read nothing from tree"));
767       continue;
768     }
769     Int_t nSdigit = fArrayS->GetEntriesFast();
770     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
771     AliInfo(Form("Got %5d digits for this event", nSdigit));
772     if (nSdigit <= 0) continue;
773     for (Int_t j = 0; j < nSdigit; j++) {
774       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
775       if (!sdigit) continue;
776       if (!ProcessSDigit(sdigit)) return kFALSE;
777     }    
778   }
779   return kTRUE;
780 }
781
782 //____________________________________________________________________
783 Bool_t 
784 AliFMDInput::ProcessRawDigits()
785 {
786   // Read the digit tree, and pass each digit to the member function
787   // ProcessDigit.
788   if (!fArrayA) {
789     AliError("No raw digit array defined");
790     return kFALSE;
791   }
792
793   Int_t nDigit = fArrayA->GetEntries();
794   if (nDigit <= 0) return kTRUE;
795   for (Int_t j = 0; j < nDigit; j++) {
796     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
797     if (!digit) continue;
798     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30) 
799       digit->Print();
800     if (!ProcessRawDigit(digit)) return kFALSE;
801   }    
802   return kTRUE;
803 }
804
805 //____________________________________________________________________
806 Bool_t 
807 AliFMDInput::ProcessRawCalibDigits()
808 {
809   // Read the digit tree, and pass each digit to the member function
810   // ProcessDigit.
811   if (!fArrayA) {
812     AliError("No raw digit array defined");
813     return kFALSE;
814   }
815
816   Int_t nDigit = fArrayA->GetEntries();
817   if (nDigit <= 0) return kTRUE;
818   for (Int_t j = 0; j < nDigit; j++) {
819     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
820     if (!digit) continue;
821     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30) 
822       digit->Print();
823     if (!ProcessRawCalibDigit(digit)) return kFALSE;
824   }    
825   return kTRUE;
826 }
827
828 //____________________________________________________________________
829 Bool_t 
830 AliFMDInput::ProcessRecPoints()
831 {
832   // Read the reconstrcted points tree, and pass each reconstruction
833   // object (AliFMDRecPoint) to either ProcessRecPoint.
834   if (!fTreeR) {
835     AliError("No recpoint tree defined");
836     return kFALSE;
837   }
838   if (!fArrayR) {
839     AliError("No recpoints array defined");
840     return kFALSE;
841   }
842
843   Int_t nEv = fTreeR->GetEntries();
844   for (Int_t i = 0; i < nEv; i++) {
845     Int_t recRead  = fTreeR->GetEntry(i);
846     if (recRead <= 0) continue;
847     Int_t nRecPoint = fArrayR->GetEntries();
848     for (Int_t j = 0; j < nRecPoint; j++) {
849       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
850       if (!recPoint) continue;
851       if (!ProcessRecPoint(recPoint)) return kFALSE;
852     }    
853   }
854   return kTRUE;
855 }
856
857 //____________________________________________________________________
858 Bool_t 
859 AliFMDInput::ProcessESDs()
860 {
861   // Process event summary data
862   if (!fESD) return kFALSE;
863   for (UShort_t det = 1; det <= 3; det++) {
864     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
865     for (Char_t* rng = rings; *rng != '\0'; rng++) {
866       UShort_t nsec = (*rng == 'I' ?  20 :  40);
867       UShort_t nstr = (*rng == 'I' ? 512 : 256);
868       for (UShort_t sec = 0; sec < nsec; sec++) {
869         for (UShort_t str = 0; str < nstr; str++) {
870           Float_t eta  = fESD->Eta(det,*rng,sec,str);
871           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
872           if (!fESD->IsAngleCorrected()) 
873             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
874           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
875         }
876       }
877     }
878   }
879   return kTRUE;
880 }
881
882 //____________________________________________________________________
883 Bool_t 
884 AliFMDInput::ProcessUsers()
885 {
886   // Process event summary data
887   for (UShort_t det = 1; det <= 3; det++) {
888     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
889     for (Char_t* rng = rings; *rng != '\0'; rng++) {
890       UShort_t nsec = (*rng == 'I' ?  20 :  40);
891       UShort_t nstr = (*rng == 'I' ? 512 : 256);
892       for (UShort_t sec = 0; sec < nsec; sec++) {
893         for (UShort_t str = 0; str < nstr; str++) {
894           Float_t v  = GetSignal(det,*rng,sec,str);
895           if (!ProcessUser(det, *rng, sec, str, v)) continue;
896         }
897       }
898     }
899   }
900   return kTRUE;
901 }
902
903 //____________________________________________________________________
904 Bool_t
905 AliFMDInput::End()
906 {
907   // Called at the end of each event.  Per default, it unloads the
908   // data trees and resets the pointers to the output arrays.   Users
909   // can overload this, but should call this member function in the
910   // overloaded member function of the derived class. 
911
912   // Check if we have been initialized
913   if (!fIsInit) { 
914     AliError("Not initialized");
915     return fIsInit;
916   }
917   // Possibly unload global kinematics information 
918   if (IsLoaded(kKinematics)) {
919     fLoader->UnloadKinematics();
920     // fTreeK = 0;
921     fStack = 0;
922   }
923   // Possibly unload FMD Hit information 
924   if (IsLoaded(kHits)) {
925     fFMDLoader->UnloadHits();
926     fTreeH = 0;
927   }
928   // Possibly unload FMD Digit information 
929   if (IsLoaded(kDigits)) {
930     fFMDLoader->UnloadDigits();
931     fTreeD = 0;
932   }
933   // Possibly unload FMD Sdigit information 
934   if (IsLoaded(kSDigits)) {
935     fFMDLoader->UnloadSDigits();
936     fTreeS = 0;
937   }
938   // Possibly unload FMD RecPoints information 
939   if (IsLoaded(kRecPoints)) {
940     fFMDLoader->UnloadRecPoints();
941     fTreeR = 0;
942   }
943   // AliInfo("Now out event");
944   return kTRUE;
945 }
946
947 //____________________________________________________________________
948 Bool_t
949 AliFMDInput::Run(UInt_t maxEvents)
950 {
951   // Run over all events and files references in galice.root 
952
953   Bool_t retval;
954   if (!(retval = Init())) return retval;
955
956   fNEvents = NEvents();
957   if (fNEvents < 0)       fNEvents = maxEvents;
958   else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
959
960   Int_t event = 0;
961   for (; fNEvents < 0 || event < fNEvents; event++) {
962     printf("\rEvent %8d/%8d ...", event, fNEvents);
963     if (!(retval = Begin(event))) break;
964     if (!(retval = Event())) break;
965     if (!(retval = End())) break;
966   }
967   printf("Looped over %8d events\n", event+1);
968   if (!retval) return retval;
969   retval = Finish();
970   return retval;
971 }
972
973 //__________________________________________________________________
974 TArrayF 
975 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max) 
976 {
977   // Service function to define a logarithmic axis. 
978   // Parameters: 
979   //   n    Number of bins 
980   //   min  Minimum of axis 
981   //   max  Maximum of axis 
982   TArrayF bins(n+1);
983   bins[0]      = min;
984   if (n <= 20) {
985     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
986     return bins;
987   }
988   Float_t dp   = n / TMath::Log10(max / min);
989   Float_t pmin = TMath::Log10(min);
990   for (Int_t i = 1; i < n+1; i++) {
991     Float_t p = pmin + i / dp;
992     bins[i]   = TMath::Power(10, p);
993   }
994   return bins;
995 }
996
997 //____________________________________________________________________
998 void
999 AliFMDInput::ScanDirectory(TSystemDirectory* dir, 
1000                            const TString& olddir, 
1001                            TChain* chain, 
1002                            const char* pattern, bool recursive)
1003 {
1004   // Get list of files, and go back to old working directory
1005   TString oldDir(gSystem->WorkingDirectory());
1006   TList* files = dir->GetListOfFiles();
1007   gSystem->ChangeDirectory(oldDir);
1008   
1009   // Sort list of files and check if we should add it 
1010   if (!files) return;
1011   files->Sort();
1012   TIter next(files);
1013   TSystemFile* file = 0;
1014   while ((file = static_cast<TSystemFile*>(next()))) {
1015     TString name(file->GetName());
1016     
1017     // Ignore special links 
1018     if (name == "." || name == "..") continue;
1019
1020     // Check if this is a directory 
1021     if (file->IsDirectory()) { 
1022       if (recursive) 
1023         ScanDirectory(static_cast<TSystemDirectory*>(file),
1024                       olddir, chain,
1025                       pattern,recursive);
1026       continue;
1027     }
1028     
1029     // If this is not a root file, ignore 
1030     if (!name.EndsWith(".root")) continue;
1031
1032     // If this file does not contain the pattern, ignore 
1033     if (!name.Contains(pattern)) continue;
1034     if (name.Contains("friends")) continue;
1035     
1036     // Get the path 
1037     TString data(Form("%s/%s", file->GetTitle(), name.Data()));
1038
1039     TFile* test = TFile::Open(data.Data(), "READ");
1040     if (!test || test->IsZombie()) { 
1041       ::Warning("ScanDirectory", "Failed to open file %s", data.Data());
1042       continue;
1043     }
1044     test->Close();
1045     chain->Add(data);
1046   }
1047 }
1048
1049 //____________________________________________________________________
1050 TChain*
1051 AliFMDInput::MakeChain(const char* what, const char* datadir, bool recursive)
1052 {
1053   TString w(what);
1054   w.ToUpper();
1055   const char* treeName = 0;
1056   const char* pattern  = 0;
1057   if      (w.Contains("ESD")) { treeName = "esdTree"; pattern = "AliESD"; }
1058   else if (w.Contains("MC"))  { treeName = "TE";      pattern = "galice"; }
1059   else {
1060     ::Error("MakeChain", "Unknown mode '%s' (not one of ESD, or MC)", what);
1061     return 0;
1062   }
1063     
1064   // --- Our data chain ----------------------------------------------
1065   TChain* chain = new TChain(treeName);
1066
1067   // --- Get list of ESDs --------------------------------------------
1068   // Open source directory, and make sure we go back to were we were 
1069   TString oldDir(gSystem->WorkingDirectory());
1070   TSystemDirectory d(datadir, datadir);
1071   ScanDirectory(&d, oldDir, chain, pattern, recursive);
1072
1073   return chain;
1074 }
1075
1076
1077 //____________________________________________________________________
1078 //
1079 // EOF
1080 //