A lot of changes after detector review:
[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 "AliFMD.h"             // ALIFMD_H
41 #include "AliFMDHit.h"          // ALIFMDHIT_H
42 #include "AliFMDDigit.h"        // ALIFMDDigit_H
43 #include "AliFMDSDigit.h"       // ALIFMDDigit_H
44 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
45 #include "AliFMDRawReader.h"    // ALIFMDRAWREADER_H
46 #include <AliESD.h>
47 #include <AliESDFMD.h>
48 #include <AliESDEvent.h>
49 #include <AliCDBManager.h>
50 #include <AliCDBEntry.h>
51 #include <AliAlignObjParams.h>
52 #include <AliTrackReference.h>
53 #include <TTree.h>              // ROOT_TTree
54 #include <TChain.h>             // ROOT_TChain
55 #include <TParticle.h>          // ROOT_TParticle
56 #include <TString.h>            // ROOT_TString
57 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
58 #include <TMath.h>              // ROOT_TMath
59 #include <TGeoManager.h>        // ROOT_TGeoManager 
60 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
61 #include <Riostream.h>          // ROOT_Riostream
62 #include <TFile.h>              // ROOT_TFile
63 #include <TStreamerInfo.h>
64 #include <TArrayF.h>
65
66 //____________________________________________________________________
67 ClassImp(AliFMDInput)
68 #if 0
69   ; // This is here to keep Emacs for indenting the next line
70 #endif
71
72
73 //____________________________________________________________________
74 AliFMDInput::AliFMDInput()
75   : TNamed("AliFMDInput", "Input handler for various FMD data"), 
76     fGAliceFile(""), 
77     fLoader(0),
78     fRun(0), 
79     fStack(0),
80     fFMDLoader(0), 
81     fReader(0),
82     fFMDReader(0),
83     fFMD(0),
84     fESD(0),
85     fESDEvent(0),
86     fTreeE(0),
87     fTreeH(0),
88     fTreeTR(0),
89     fTreeD(0),
90     fTreeS(0),
91     fTreeR(0), 
92     fTreeA(0), 
93     fChainE(0),
94     fArrayE(0),
95     fArrayH(0),
96     fArrayTR(0), 
97     fArrayD(0),
98     fArrayS(0), 
99     fArrayR(0), 
100     fArrayA(0), 
101     fHeader(0),
102     fGeoManager(0),
103     fTreeMask(0), 
104     fRawFile(""),
105     fIsInit(kFALSE),
106     fEventCount(0)
107 {
108
109   // Constructor of an FMD input object.  Specify what data to read in
110   // using the AddLoad member function.   Sub-classes should at a
111   // minimum overload the member function Event.   A full job can be
112   // executed using the member function Run. 
113 }
114
115   
116
117 //____________________________________________________________________
118 AliFMDInput::AliFMDInput(const char* gAliceFile)
119   : TNamed("AliFMDInput", "Input handler for various FMD data"), 
120     fGAliceFile(gAliceFile),
121     fLoader(0),
122     fRun(0), 
123     fStack(0),
124     fFMDLoader(0), 
125     fReader(0),
126     fFMDReader(0),
127     fFMD(0),
128     fESD(0),
129     fESDEvent(0),
130     fTreeE(0),
131     fTreeH(0),
132     fTreeTR(0),
133     fTreeD(0),
134     fTreeS(0),
135     fTreeR(0), 
136     fTreeA(0), 
137     fChainE(0),
138     fArrayE(0),
139     fArrayH(0),
140     fArrayTR(0), 
141     fArrayD(0),
142     fArrayS(0), 
143     fArrayR(0), 
144     fArrayA(0), 
145     fHeader(0),
146     fGeoManager(0),
147     fTreeMask(0), 
148     fRawFile(""),
149     fIsInit(kFALSE),
150     fEventCount(0)
151 {
152   
153   // Constructor of an FMD input object.  Specify what data to read in
154   // using the AddLoad member function.   Sub-classes should at a
155   // minimum overload the member function Event.   A full job can be
156   // executed using the member function Run. 
157 }
158
159 //____________________________________________________________________
160 Int_t
161 AliFMDInput::NEvents() const 
162 {
163   // Get number of events
164   if (TESTBIT(fTreeMask, kRaw)) return fReader->GetNumberOfEvents();
165   if (fChainE) return fChainE->GetEntriesFast();
166   if (fTreeE) return fTreeE->GetEntries();
167   return -1;
168 }
169
170 //____________________________________________________________________
171 Bool_t
172 AliFMDInput::Init()
173 {
174   // Initialize the object.  Get the needed loaders, and such. 
175
176   // Check if we have been initialized
177   if (fIsInit) { 
178     AliWarning("Already initialized");
179     return fIsInit;
180   }
181   Info("Init","Initialising w/mask 0x%04x\n"
182        "\tHits:          %d\n"
183        "\tKinematics:    %d\n"
184        "\tDigits:        %d\n"
185        "\tSDigits:       %d\n"
186        "\tHeader:        %d\n"
187        "\tRecPoints:     %d\n"
188        "\tESD:           %d\n"
189        "\tRaw:           %d\n"
190        "\tGeometry:      %d\n"
191        "\tTracks:        %d\n"
192        "\tTracksRefs:    %d",
193        fTreeMask,
194        TESTBIT(fTreeMask, kHits),
195        TESTBIT(fTreeMask, kKinematics),
196        TESTBIT(fTreeMask, kDigits),
197        TESTBIT(fTreeMask, kSDigits),
198        TESTBIT(fTreeMask, kHeader),
199        TESTBIT(fTreeMask, kRecPoints),
200        TESTBIT(fTreeMask, kESD),
201        TESTBIT(fTreeMask, kRaw),
202        TESTBIT(fTreeMask, kGeometry),
203        TESTBIT(fTreeMask, kTracks),
204        TESTBIT(fTreeMask, kTrackRefs));
205   // Get the run 
206   if (TESTBIT(fTreeMask, kDigits)     ||
207       TESTBIT(fTreeMask, kSDigits)    || 
208       TESTBIT(fTreeMask, kKinematics) || 
209       TESTBIT(fTreeMask, kTrackRefs)  || 
210       TESTBIT(fTreeMask, kTracks)     || 
211       TESTBIT(fTreeMask, kHeader)) {
212     if (!gSystem->FindFile(".:/", fGAliceFile)) {
213       AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
214     }
215     else {
216       fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
217       if (!fLoader) {
218         AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
219         return kFALSE;
220       }
221       AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
222
223       if  (fLoader->LoadgAlice()) return kFALSE;
224       
225       fRun = fLoader->GetAliRun();
226       
227       // Get the FMD 
228       fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
229       if (!fFMD) {
230         AliError("Failed to get detector FMD from loader");
231         return kFALSE;
232       }
233       
234       // Get the FMD loader
235       fFMDLoader = fLoader->GetLoader("FMDLoader");
236       if (!fFMDLoader) {
237         AliError("Failed to get detector FMD loader from loader");
238         return kFALSE;
239       }
240       if (fLoader->LoadHeader()) { 
241         AliError("Failed to get event header information from loader");
242         return kFALSE;
243       }
244       fTreeE = fLoader->TreeE();
245     }
246   }
247
248   // Optionally, get the ESD files
249   if (TESTBIT(fTreeMask, kESD)) {
250     fChainE = new TChain("esdTree");
251     TSystemDirectory dir(".",".");
252     TList*           files = dir.GetListOfFiles();
253     TSystemFile*     file = 0;
254     if (!files) {
255       AliError("No files");
256       return kFALSE;
257     }
258     files->Sort();
259     TIter            next(files);
260     while ((file = static_cast<TSystemFile*>(next()))) {
261       TString fname(file->GetName());
262       if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
263     }
264     fESDEvent = new AliESDEvent();
265     fESDEvent->ReadFromTree(fChainE);
266     //    fChainE->SetBranchAddress("ESD", &fMainESD);
267     
268   }
269     
270   if (TESTBIT(fTreeMask, kRaw)) {
271     AliInfo("Getting FMD raw data digits");
272     fArrayA = new TClonesArray("AliFMDDigit");
273 #if 0
274     if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
275       fReader = new AliRawReaderRoot(fRawFile.Data());
276     else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
277       fReader = new AliRawReaderDate(fRawFile.Data());
278     else
279       fReader = new AliRawReaderFile(-1);
280 #else
281     if(!fRawFile.IsNull()) 
282       fReader = AliRawReader::Create(fRawFile.Data());
283     else 
284       fReader = new AliRawReaderFile(-1);
285 #endif
286     fFMDReader = new AliFMDRawReader(fReader, 0);
287   }
288   
289   // Optionally, get the geometry 
290   if (TESTBIT(fTreeMask, kGeometry)) {
291     if (fRun) {
292       TString fname(fRun->GetGeometryFileName());
293       if (fname.IsNull()) {
294         Warning("Init", "No file name for the geometry from AliRun");
295         fname = gSystem->DirName(fGAliceFile);
296         fname.Append("/geometry.root");
297       }
298       fGeoManager = TGeoManager::Import(fname.Data());
299       if (!fGeoManager) {
300         Fatal("Init", "No geometry manager found");
301         return kFALSE;
302       }
303     }
304     else { 
305       AliGeomManager::LoadGeometry();
306     }
307     AliCDBManager* cdb   = AliCDBManager::Instance();
308     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
309     if (align) {
310       AliInfo("Got alignment data from CDB");
311       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
312       if (!array) {
313         AliWarning("Invalid align data from CDB");
314       }
315       else {
316         Int_t nAlign = array->GetEntries();
317         for (Int_t i = 0; i < nAlign; i++) {
318           AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
319           if (!a->ApplyToGeometry()) {
320             AliWarning(Form("Failed to apply alignment to %s", 
321                             a->GetSymName()));
322           }
323         }
324       }
325     }
326   }
327
328   fEventCount = 0;
329   fIsInit = kTRUE;
330   return fIsInit;
331 }
332
333 //____________________________________________________________________
334 Bool_t
335 AliFMDInput::Begin(Int_t event)
336 {
337   // Called at the begining of each event.  Per default, it gets the
338   // data trees and gets pointers to the output arrays.   Users can
339   // overload this, but should call this member function in the
340   // overloaded member function of the derived class. 
341
342   // Check if we have been initialized
343   if (!fIsInit) { 
344     AliError("Not initialized");
345     return fIsInit;
346   }
347
348   // Get the event 
349   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
350   AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
351
352   // Possibly load global kinematics information 
353   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
354     // AliInfo("Getting kinematics");
355     if (fLoader->LoadKinematics("READ")) return kFALSE;
356     fStack = fLoader->Stack();
357   }
358
359   // Possibly load FMD Hit information 
360   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
361     // AliInfo("Getting FMD hits");
362     if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
363     fTreeH = fFMDLoader->TreeH();
364     if (!fArrayH) fArrayH = fFMD->Hits(); 
365   }
366   
367   // Possibly load FMD TrackReference information 
368   if (TESTBIT(fTreeMask, kTrackRefs) || TESTBIT(fTreeMask, kTracks)) {
369     // AliInfo("Getting FMD hits");
370     if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
371     fTreeTR = fLoader->TreeTR();
372     if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
373     fTreeTR->SetBranchAddress("TrackReferences",  &fArrayTR);
374   }
375   
376   // Possibly load heaedr information 
377   if (TESTBIT(fTreeMask, kHeader)) {
378     // AliInfo("Getting FMD hits");
379     if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
380     fHeader = fLoader->GetHeader();
381   }
382
383   // Possibly load FMD Digit information 
384   if (TESTBIT(fTreeMask, kDigits)) {
385     // AliInfo("Getting FMD digits");
386     if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
387     fTreeD = fFMDLoader->TreeD();
388     if (fTreeD) {
389       if (!fArrayD) fArrayD = fFMD->Digits();
390     }
391     else {
392       fArrayD = 0;
393       AliWarning(Form("Failed to load FMD Digits"));
394     } 
395   }
396
397   // Possibly load FMD Sdigit information 
398   if (TESTBIT(fTreeMask, kSDigits)) {
399     // AliInfo("Getting FMD summable digits");
400     if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) { 
401       AliWarning("Failed to load SDigits!");
402       return kFALSE;
403     }
404     fTreeS = fFMDLoader->TreeS();
405     if (!fArrayS) fArrayS = fFMD->SDigits();
406   }
407
408   // Possibly load FMD RecPoints information 
409   if (TESTBIT(fTreeMask, kRecPoints)) {
410     // AliInfo("Getting FMD reconstructed points");
411     if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
412     fTreeR = fFMDLoader->TreeR();
413     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
414     fTreeR->SetBranchAddress("FMD",  &fArrayR);
415   }  
416
417   // Possibly load FMD ESD information 
418   if (TESTBIT(fTreeMask, kESD)) {
419     // AliInfo("Getting FMD event summary data");
420     Int_t read = fChainE->GetEntry(event);
421     if (read <= 0) return kFALSE;
422     fESD = fESDEvent->GetFMDData();
423     if (!fESD) return kFALSE;
424 #if 0
425     TFile* f = fChainE->GetFile();
426     if (f) {
427       TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
428       if (o) {
429         TStreamerInfo* info = static_cast<TStreamerInfo*>(o);
430         std::cout << "AliFMDMap class version read is " 
431                   <<  info->GetClassVersion() << std::endl;
432       }
433     }
434     // fESD->CheckNeedUShort(fChainE->GetFile());
435 #endif
436   }
437
438   // Possibly load FMD Digit information 
439   if (TESTBIT(fTreeMask, kRaw)) {
440     // AliInfo("Getting FMD raw data digits");
441     if (!fReader->NextEvent()) return kFALSE;
442     // AliFMDRawReader r(fReader, 0);
443     fArrayA->Clear();
444     fFMDReader->ReadAdcs(fArrayA);
445     AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
446   }
447   fEventCount++;
448   return kTRUE;
449 }
450
451
452 //____________________________________________________________________
453 Bool_t 
454 AliFMDInput::Event()
455 {
456   // Process one event.  The default implementation one or more of 
457   //
458   //   -  ProcessHits       if the hits are loaded. 
459   //   -  ProcessDigits     if the digits are loaded. 
460   //   -  ProcessSDigits    if the sumbable digits are loaded. 
461   //   -  ProcessRecPoints  if the reconstructed points are loaded. 
462   //   -  ProcessESD        if the event summary data is loaded
463   // 
464   if (TESTBIT(fTreeMask, kHits)) 
465     if (!ProcessHits()) return kFALSE; 
466   if (TESTBIT(fTreeMask, kTrackRefs)) 
467     if (!ProcessTrackRefs()) return kFALSE; 
468   if (TESTBIT(fTreeMask, kTracks)) 
469     if (!ProcessTracks()) return kFALSE; 
470   if (TESTBIT(fTreeMask, kSDigits)) 
471     if (!ProcessSDigits()) return kFALSE;
472   if (TESTBIT(fTreeMask, kDigits)) 
473     if (!ProcessDigits()) return kFALSE;
474   if (TESTBIT(fTreeMask, kRaw)) 
475     if (!ProcessRawDigits()) return kFALSE;
476   if (TESTBIT(fTreeMask, kRecPoints)) 
477     if (!ProcessRecPoints()) return kFALSE;
478   if (TESTBIT(fTreeMask, kESD))
479     if (!ProcessESDs()) return kFALSE;
480   
481   return kTRUE;
482 }
483
484 //____________________________________________________________________
485 Bool_t 
486 AliFMDInput::ProcessHits()
487 {
488   // Read the hit tree, and pass each hit to the member function
489   // ProcessHit.
490   if (!fTreeH) {
491     AliError("No hit tree defined");
492     return kFALSE;
493   }
494   if (!fArrayH) {
495     AliError("No hit array defined");
496     return kFALSE;
497   }
498
499   Int_t nTracks = fTreeH->GetEntries();
500   for (Int_t i = 0; i < nTracks; i++) {
501     Int_t hitRead  = fTreeH->GetEntry(i);
502     if (hitRead <= 0) continue;
503
504     Int_t nHit = fArrayH->GetEntries();
505     if (nHit <= 0) continue;
506   
507     for (Int_t j = 0; j < nHit; j++) {
508       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
509       if (!hit) continue;
510
511       TParticle* track = 0;
512       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
513         Int_t trackno = hit->Track();
514         track = fStack->Particle(trackno);
515       }
516       if (!ProcessHit(hit, track)) return kFALSE;
517     }    
518   }
519   return kTRUE;
520 }
521 //____________________________________________________________________
522 Bool_t 
523 AliFMDInput::ProcessTrackRefs()
524 {
525   // Read the reconstrcted points tree, and pass each reconstruction
526   // object (AliFMDRecPoint) to either ProcessRecPoint.
527   if (!fTreeTR) {
528     AliError("No track reference tree defined");
529     return kFALSE;
530   }
531   if (!fArrayTR) {
532     AliError("No track reference array defined");
533     return kFALSE;
534   }
535
536   Int_t nEv = fTreeTR->GetEntries();
537   for (Int_t i = 0; i < nEv; i++) {
538     Int_t trRead  = fTreeTR->GetEntry(i);
539     if (trRead <= 0) continue;
540     Int_t nTrackRefs = fArrayTR->GetEntries();
541     for (Int_t j = 0; j < nTrackRefs; j++) {
542       AliTrackReference* trackRef = static_cast<AliTrackReference*>(fArrayTR->At(j));
543       if (!trackRef) continue;
544       if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
545       TParticle* track = 0;
546       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
547         Int_t trackno = trackRef->GetTrack();
548         track = fStack->Particle(trackno);
549       }
550       if (!ProcessTrackRef(trackRef,track)) return kFALSE;
551     }    
552   }
553   return kTRUE;
554 }
555 //____________________________________________________________________
556 Bool_t 
557 AliFMDInput::ProcessTracks()
558 {
559   // Read the hit tree, and pass each hit to the member function
560   // ProcessTrack.
561   if (!fStack) {
562     AliError("No track tree defined");
563     return kFALSE;
564   }
565   if (!fTreeH) {
566     AliError("No hit tree defined");
567     return kFALSE;
568   }
569   if (!fArrayH) {
570     AliError("No hit array defined");
571     return kFALSE;
572   }
573
574   // Int_t nTracks = fStack->GetNtrack();
575   Int_t nTracks = fTreeH->GetEntries();
576   for (Int_t i = 0; i < nTracks; i++) {
577     Int_t      trackno = nTracks - i - 1;
578     TParticle* track   = fStack->Particle(trackno);
579     if (!track) continue;
580
581     // Get the hits for this track. 
582     Int_t hitRead  = fTreeH->GetEntry(i);
583     Int_t nHit     = fArrayH->GetEntries();
584     if (nHit == 0 || hitRead <= 0) { 
585       // Let user code see the track, even if there's no hits. 
586       if (!ProcessTrack(trackno, track, 0)) return kFALSE;
587       continue;
588     }
589
590     // Loop over the hits corresponding to this track.
591     for (Int_t j = 0; j < nHit; j++) {
592       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
593       if (!ProcessTrack(trackno, track, hit)) return kFALSE;
594     }   
595   }
596   return kTRUE;
597 }
598
599 //____________________________________________________________________
600 Bool_t 
601 AliFMDInput::ProcessDigits()
602 {
603   // Read the digit tree, and pass each digit to the member function
604   // ProcessDigit.
605   if (!fTreeD) {
606     AliError("No digit tree defined");
607     return kFALSE;
608   }
609   if (!fArrayD) {
610     AliError("No digit array defined");
611     return kFALSE;
612   }
613
614   Int_t nEv = fTreeD->GetEntries();
615   for (Int_t i = 0; i < nEv; i++) {
616     Int_t digitRead  = fTreeD->GetEntry(i);
617     if (digitRead <= 0) continue;
618     Int_t nDigit = fArrayD->GetEntries();
619     AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
620     if (nDigit <= 0) continue;
621     for (Int_t j = 0; j < nDigit; j++) {
622       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
623       if (!digit) continue;
624       if (!ProcessDigit(digit)) return kFALSE;
625     }    
626   }
627   return kTRUE;
628 }
629
630 //____________________________________________________________________
631 Bool_t 
632 AliFMDInput::ProcessSDigits()
633 {
634   // Read the summable digit tree, and pass each sumable digit to the
635   // member function ProcessSdigit.
636   if (!fTreeS) {
637     AliWarning("No sdigit tree defined");
638     return kTRUE; // Empty SDigits is fine
639   }
640   if (!fArrayS) {
641     AliWarning("No sdigit array defined");
642     return kTRUE; // Empty SDigits is fine
643   }
644
645   Int_t nEv = fTreeS->GetEntries();
646   for (Int_t i = 0; i < nEv; i++) {
647     Int_t sdigitRead  = fTreeS->GetEntry(i);
648     if (sdigitRead <= 0) { 
649       AliInfo(Form("Read nothing from tree"));
650       continue;
651     }
652     Int_t nSdigit = fArrayS->GetEntriesFast();
653     AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
654     AliInfo(Form("Got %5d digits for this event", nSdigit));
655     if (nSdigit <= 0) continue;
656     for (Int_t j = 0; j < nSdigit; j++) {
657       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
658       if (!sdigit) continue;
659       if (!ProcessSDigit(sdigit)) return kFALSE;
660     }    
661   }
662   return kTRUE;
663 }
664
665 //____________________________________________________________________
666 Bool_t 
667 AliFMDInput::ProcessRawDigits()
668 {
669   // Read the digit tree, and pass each digit to the member function
670   // ProcessDigit.
671   if (!fArrayA) {
672     AliError("No raw digit array defined");
673     return kFALSE;
674   }
675
676   Int_t nDigit = fArrayA->GetEntries();
677   if (nDigit <= 0) return kTRUE;
678   for (Int_t j = 0; j < nDigit; j++) {
679     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
680     if (!digit) continue;
681     if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30) 
682       digit->Print();
683     if (!ProcessRawDigit(digit)) return kFALSE;
684   }    
685   return kTRUE;
686 }
687
688 //____________________________________________________________________
689 Bool_t 
690 AliFMDInput::ProcessRecPoints()
691 {
692   // Read the reconstrcted points tree, and pass each reconstruction
693   // object (AliFMDRecPoint) to either ProcessRecPoint.
694   if (!fTreeR) {
695     AliError("No recpoint tree defined");
696     return kFALSE;
697   }
698   if (!fArrayR) {
699     AliError("No recpoints array defined");
700     return kFALSE;
701   }
702
703   Int_t nEv = fTreeR->GetEntries();
704   for (Int_t i = 0; i < nEv; i++) {
705     Int_t recRead  = fTreeR->GetEntry(i);
706     if (recRead <= 0) continue;
707     Int_t nRecPoint = fArrayR->GetEntries();
708     for (Int_t j = 0; j < nRecPoint; j++) {
709       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
710       if (!recPoint) continue;
711       if (!ProcessRecPoint(recPoint)) return kFALSE;
712     }    
713   }
714   return kTRUE;
715 }
716
717 //____________________________________________________________________
718 Bool_t 
719 AliFMDInput::ProcessESDs()
720 {
721   // Process event summary data
722   if (!fESD) return kFALSE;
723   for (UShort_t det = 1; det <= 3; det++) {
724     Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
725     for (Char_t* rng = rings; *rng != '\0'; rng++) {
726       UShort_t nsec = (*rng == 'I' ?  20 :  40);
727       UShort_t nstr = (*rng == 'I' ? 512 : 256);
728       for (UShort_t sec = 0; sec < nsec; sec++) {
729         for (UShort_t str = 0; str < nstr; str++) {
730           Float_t eta  = fESD->Eta(det,*rng,sec,str);
731           Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
732           if (!fESD->IsAngleCorrected()) 
733             mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
734           if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
735         }
736       }
737     }
738   }
739   return kTRUE;
740 }
741
742 //____________________________________________________________________
743 Bool_t
744 AliFMDInput::End()
745 {
746   // Called at the end of each event.  Per default, it unloads the
747   // data trees and resets the pointers to the output arrays.   Users
748   // can overload this, but should call this member function in the
749   // overloaded member function of the derived class. 
750
751   // Check if we have been initialized
752   if (!fIsInit) { 
753     AliError("Not initialized");
754     return fIsInit;
755   }
756   // Possibly unload global kinematics information 
757   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
758     fLoader->UnloadKinematics();
759     // fTreeK = 0;
760     fStack = 0;
761   }
762   // Possibly unload FMD Hit information 
763   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
764     fFMDLoader->UnloadHits();
765     fTreeH = 0;
766   }
767   // Possibly unload FMD Digit information 
768   if (TESTBIT(fTreeMask, kDigits)) {
769     fFMDLoader->UnloadDigits();
770     fTreeD = 0;
771   }
772   // Possibly unload FMD Sdigit information 
773   if (TESTBIT(fTreeMask, kSDigits)) {
774     fFMDLoader->UnloadSDigits();
775     fTreeS = 0;
776   }
777   // Possibly unload FMD RecPoints information 
778   if (TESTBIT(fTreeMask, kRecPoints)) {
779     fFMDLoader->UnloadRecPoints();
780     fTreeR = 0;
781   }
782   // AliInfo("Now out event");
783   return kTRUE;
784 }
785
786 //____________________________________________________________________
787 Bool_t
788 AliFMDInput::Run()
789 {
790   // Run over all events and files references in galice.root 
791
792   Bool_t retval;
793   if (!(retval = Init())) return retval;
794
795   Int_t nEvents = NEvents();
796   for (Int_t event = 0; nEvents < 0 || event < nEvents; event++) {
797     if (!(retval = Begin(event))) break;
798     if (!(retval = Event())) break;
799     if (!(retval = End())) break;
800   }
801   if (!retval) return retval;
802   retval = Finish();
803   return retval;
804 }
805
806 //__________________________________________________________________
807 TArrayF 
808 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max) 
809 {
810   // Service function to define a logarithmic axis. 
811   // Parameters: 
812   //   n    Number of bins 
813   //   min  Minimum of axis 
814   //   max  Maximum of axis 
815   TArrayF bins(n+1);
816   bins[0]      = min;
817   if (n <= 20) {
818     for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
819     return bins;
820   }
821   Float_t dp   = n / TMath::Log10(max / min);
822   Float_t pmin = TMath::Log10(min);
823   for (Int_t i = 1; i < n+1; i++) {
824     Float_t p = pmin + i / dp;
825     bins[i]   = TMath::Power(10, p);
826   }
827   return bins;
828 }
829
830
831
832 //____________________________________________________________________
833 //
834 // EOF
835 //