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