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