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