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