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