]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - FMD/AliFMDInput.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Added new 'status' histogram to the DQM data maker and checker. This histogram
[u/mrichter/AliRoot.git] / FMD / AliFMDInput.cxx
... / ...
CommitLineData
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//____________________________________________________________________
69ClassImp(AliFMDInput)
70#if 0
71 ; // This is here to keep Emacs for indenting the next line
72#endif
73
74//____________________________________________________________________
75const AliFMDInput::ETrees AliFMDInput::fgkAllLoads[] = { kHits,
76 kKinematics,
77 kDigits,
78 kSDigits,
79 kHeader,
80 kRecPoints,
81 kESD,
82 kRaw,
83 kGeometry,
84 kTrackRefs,
85 kRawCalib,
86 kUser };
87
88//____________________________________________________________________
89AliFMDInput::AliFMDInput()
90 : TNamed("AliFMDInput", "Input handler for various FMD data"),
91 fGAliceFile(""),
92 fLoader(0),
93 fRun(0),
94 fStack(0),
95 fFMDLoader(0),
96 fReader(0),
97 fFMDReader(0),
98 fFMD(0),
99 fESD(0),
100 fESDEvent(0),
101 fTreeE(0),
102 fTreeH(0),
103 fTreeTR(0),
104 fTreeD(0),
105 fTreeS(0),
106 fTreeR(0),
107 fTreeA(0),
108 fChainE(0),
109 fArrayE(0),
110 fArrayH(0),
111 fArrayTR(0),
112 fArrayD(0),
113 fArrayS(0),
114 fArrayR(0),
115 fArrayA(0),
116 fHeader(0),
117 fGeoManager(0),
118 fTreeMask(0),
119 fRawFile(""),
120 fInputDir("."),
121 fIsInit(kFALSE),
122 fEventCount(0),
123 fNEvents(-1)
124{
125
126 // Constructor of an FMD input object. Specify what data to read in
127 // using the AddLoad member function. Sub-classes should at a
128 // minimum overload the member function Event. A full job can be
129 // executed using the member function Run.
130}
131
132
133
134//____________________________________________________________________
135AliFMDInput::AliFMDInput(const char* gAliceFile)
136 : TNamed("AliFMDInput", "Input handler for various FMD data"),
137 fGAliceFile(gAliceFile),
138 fLoader(0),
139 fRun(0),
140 fStack(0),
141 fFMDLoader(0),
142 fReader(0),
143 fFMDReader(0),
144 fFMD(0),
145 fESD(0),
146 fESDEvent(0),
147 fTreeE(0),
148 fTreeH(0),
149 fTreeTR(0),
150 fTreeD(0),
151 fTreeS(0),
152 fTreeR(0),
153 fTreeA(0),
154 fChainE(0),
155 fArrayE(0),
156 fArrayH(0),
157 fArrayTR(0),
158 fArrayD(0),
159 fArrayS(0),
160 fArrayR(0),
161 fArrayA(0),
162 fHeader(0),
163 fGeoManager(0),
164 fTreeMask(0),
165 fRawFile(""),
166 fInputDir("."),
167 fIsInit(kFALSE),
168 fEventCount(0),
169 fNEvents(-1)
170{
171
172 // Constructor of an FMD input object. Specify what data to read in
173 // using the AddLoad member function. Sub-classes should at a
174 // minimum overload the member function Event. A full job can be
175 // executed using the member function Run.
176}
177
178//____________________________________________________________________
179void
180AliFMDInput::SetLoads(UInt_t mask)
181{
182 for (UInt_t i = 0; i < sizeof(mask); i++) {
183 if (!(mask & (1 << i))) continue;
184 const ETrees *ptype = fgkAllLoads;
185 do {
186 ETrees type = *ptype;
187 if (i != UInt_t(type)) continue;
188 AddLoad(type);
189 break;
190 } while (*ptype++ != kUser);
191 }
192}
193
194//____________________________________________________________________
195void
196AliFMDInput::SetLoads(const char* what)
197{
198 TString l(what);
199 TObjArray* ll = l.Tokenize(", ");
200 TIter next(ll);
201 TObject* os = 0;
202 while ((os = next())) {
203 ETrees type = ParseLoad(os->GetName());
204 AddLoad(type);
205 }
206}
207
208
209//____________________________________________________________________
210AliFMDInput::ETrees
211AliFMDInput::ParseLoad(const char* what)
212{
213 TString opt(what);
214 opt.ToLower();
215 const ETrees* ptype = fgkAllLoads;
216 do {
217 ETrees type = *ptype;
218 if (opt.Contains(TreeName(type,true), TString::kIgnoreCase))
219 return type;
220 } while (*ptype++ != kUser);
221 return kUser;
222}
223//____________________________________________________________________
224const char*
225AliFMDInput::LoadedString(Bool_t dataOnly) const
226{
227 static TString ret;
228 if (!ret.IsNull()) return ret.Data();
229
230 const ETrees* ptype = fgkAllLoads;
231 do {
232 ETrees type = *ptype;
233 if (dataOnly &&
234 (type == kKinematics ||
235 type == kHeader ||
236 type == kGeometry ||
237 type == kTrackRefs)) continue;
238 if (!IsLoaded(*ptype)) continue;
239
240 if (!ret.IsNull()) ret.Append(",");
241 ret.Append(TreeName(type));
242 } while (*ptype++ != kUser);
243 return ret.Data();
244}
245
246//____________________________________________________________________
247const char*
248AliFMDInput::TreeName(ETrees tree, Bool_t shortest)
249{
250 if (shortest) {
251 switch (tree) {
252 case kHits: return "hit";
253 case kKinematics: return "kin";
254 case kDigits: return "dig";
255 case kSDigits: return "sdig";
256 case kHeader: return "hea";
257 case kRecPoints: return "recp";
258 case kESD: return "esd";
259 case kRaw: return "raw";
260 case kGeometry: return "geo";
261 case kTrackRefs: return "trackr";
262 case kRawCalib: return "rawc";
263 case kUser: return "user";
264 }
265 return 0;
266 }
267 switch (tree) {
268 case kHits: return "Hits";
269 case kKinematics: return "Kinematics";
270 case kDigits: return "Digits";
271 case kSDigits: return "SDigits";
272 case kHeader: return "Header";
273 case kRecPoints: return "RecPoints";
274 case kESD: return "ESD";
275 case kRaw: return "Raw";
276 case kGeometry: return "Geometry";
277 case kTrackRefs: return "TrackRefs";
278 case kRawCalib: return "RawCalib";
279 case kUser: return "User";
280 }
281 return 0;
282}
283
284//____________________________________________________________________
285Int_t
286AliFMDInput::NEvents() const
287{
288 // Get number of events
289 if (IsLoaded(kRaw) ||
290 IsLoaded(kRawCalib)) return fReader->GetNumberOfEvents();
291 if (fChainE) return fChainE->GetEntriesFast();
292 if (fTreeE) return fTreeE->GetEntries();
293 return -1;
294}
295
296//____________________________________________________________________
297Bool_t
298AliFMDInput::Init()
299{
300 // Initialize the object. Get the needed loaders, and such.
301
302 // Check if we have been initialized
303 if (fIsInit) {
304 AliWarning("Already initialized");
305 return fIsInit;
306 }
307 TString what;
308 const ETrees* ptype = fgkAllLoads;
309 do {
310 ETrees type = *ptype;
311 what.Append(Form("\n\t%-20s: %s", TreeName(type),
312 IsLoaded(type) ? "yes" : "no"));
313 } while (*ptype++ != kUser);
314
315 Info("Init","Initialising w/mask 0x%04x%s", fTreeMask, what.Data());
316 // Get the run
317 if (IsLoaded(kDigits) ||
318 IsLoaded(kSDigits) ||
319 IsLoaded(kKinematics) ||
320 IsLoaded(kTrackRefs) ||
321 IsLoaded(kHeader)) {
322 if (!gSystem->FindFile(".:/", fGAliceFile)) {
323 AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
324 }
325 else {
326 fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
327 if (!fLoader) {
328 AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
329 return kFALSE;
330 }
331 AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
332
333 if (fLoader->LoadgAlice()) return kFALSE;
334
335 fRun = fLoader->GetAliRun();
336
337 // Get the FMD
338 fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
339 if (!fFMD) {
340 AliError("Failed to get detector FMD from loader");
341 return kFALSE;
342 }
343
344 // Get the FMD loader
345 fFMDLoader = fLoader->GetLoader("FMDLoader");
346 if (!fFMDLoader) {
347 AliError("Failed to get detector FMD loader from loader");
348 return kFALSE;
349 }
350 if (fLoader->LoadHeader()) {
351 AliError("Failed to get event header information from loader");
352 return kFALSE;
353 }
354 fTreeE = fLoader->TreeE();
355 }
356 }
357
358 // Optionally, get the ESD files
359 if (IsLoaded(kESD)) {
360 fChainE = MakeChain("ESD", fInputDir, true);
361 fESDEvent = new AliESDEvent();
362 fESDEvent->ReadFromTree(fChainE);
363 // fChainE->SetBranchAddress("ESD", &fMainESD);
364
365 }
366
367 if (IsLoaded(kRaw) ||
368 IsLoaded(kRawCalib)) {
369 AliInfo("Getting FMD raw data digits");
370 fArrayA = new TClonesArray("AliFMDDigit");
371#if 0
372 if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
373 fReader = new AliRawReaderRoot(fRawFile.Data());
374 else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
375 fReader = new AliRawReaderDate(fRawFile.Data());
376 else
377 fReader = new AliRawReaderFile(-1);
378#else
379 if(!fRawFile.IsNull())
380 fReader = AliRawReader::Create(fRawFile.Data());
381 else
382 fReader = new AliRawReaderFile(-1);
383#endif
384 fFMDReader = new AliFMDRawReader(fReader, 0);
385 }
386
387 // Optionally, get the geometry
388 if (IsLoaded(kGeometry)) {
389 TString fname;
390 if (fRun) {
391 fname = gSystem->DirName(fGAliceFile);
392 fname.Append("/geometry.root");
393 }
394 if (!gSystem->AccessPathName(fname.Data()))
395 fname = "";
396 AliCDBManager* cdb = AliCDBManager::Instance();
397 if (!cdb->IsDefaultStorageSet()) {
398 cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
399 cdb->SetRun(0);
400 }
401
402 AliGeomManager::LoadGeometry(fname.IsNull() ? 0 : fname.Data());
403
404 AliCDBEntry* align = cdb->Get("FMD/Align/Data");
405 if (align) {
406 AliInfo("Got alignment data from CDB");
407 TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
408 if (!array) {
409 AliWarning("Invalid align data from CDB");
410 }
411 else {
412 Int_t nAlign = array->GetEntries();
413 for (Int_t i = 0; i < nAlign; i++) {
414 AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
415 if (!a->ApplyToGeometry()) {
416 AliWarning(Form("Failed to apply alignment to %s",
417 a->GetSymName()));
418 }
419 }
420 }
421 }
422 AliFMDGeometry* geom = AliFMDGeometry::Instance();
423 geom->Init();
424 geom->InitTransformations();
425 }
426
427 fEventCount = 0;
428 fIsInit = kTRUE;
429 return fIsInit;
430}
431
432//____________________________________________________________________
433Bool_t
434AliFMDInput::Begin(Int_t event)
435{
436 // Called at the begining of each event. Per default, it gets the
437 // data trees and gets pointers to the output arrays. Users can
438 // overload this, but should call this member function in the
439 // overloaded member function of the derived class.
440
441 // Check if we have been initialized
442 if (!fIsInit) {
443 AliError("Not initialized");
444 return fIsInit;
445 }
446
447 // Get the event
448 if (fLoader && fLoader->GetEvent(event)) return kFALSE;
449
450 // Possibly load global kinematics information
451 if (IsLoaded(kKinematics)) {
452 // AliInfo("Getting kinematics");
453 if (fLoader->LoadKinematics("READ")) return kFALSE;
454 fStack = fLoader->Stack();
455 }
456
457 // Possibly load FMD Hit information
458 if (IsLoaded(kHits)) {
459 // AliInfo("Getting FMD hits");
460 if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
461 fTreeH = fFMDLoader->TreeH();
462 if (!fArrayH) fArrayH = fFMD->Hits();
463 }
464
465 // Possibly load FMD TrackReference information
466 if (IsLoaded(kTrackRefs)) {
467 // AliInfo("Getting FMD hits");
468 if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
469 fTreeTR = fLoader->TreeTR();
470 if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
471 fTreeTR->SetBranchAddress("TrackReferences", &fArrayTR);
472 }
473
474 // Possibly load heaedr information
475 if (IsLoaded(kHeader)) {
476 // AliInfo("Getting FMD hits");
477 if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
478 fHeader = fLoader->GetHeader();
479 }
480
481 // Possibly load FMD Digit information
482 if (IsLoaded(kDigits)) {
483 // AliInfo("Getting FMD digits");
484 if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
485 fTreeD = fFMDLoader->TreeD();
486 if (fTreeD) {
487 if (!fArrayD) fArrayD = fFMD->Digits();
488 }
489 else {
490 fArrayD = 0;
491 AliWarning(Form("Failed to load FMD Digits"));
492 }
493 }
494
495 // Possibly load FMD Sdigit information
496 if (IsLoaded(kSDigits)) {
497 // AliInfo("Getting FMD summable digits");
498 if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) {
499 AliWarning("Failed to load SDigits!");
500 return kFALSE;
501 }
502 fTreeS = fFMDLoader->TreeS();
503 if (!fArrayS) fArrayS = fFMD->SDigits();
504 }
505
506 // Possibly load FMD RecPoints information
507 if (IsLoaded(kRecPoints)) {
508 // AliInfo("Getting FMD reconstructed points");
509 if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
510 fTreeR = fFMDLoader->TreeR();
511 if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
512 fTreeR->SetBranchAddress("FMD", &fArrayR);
513 }
514
515 // Possibly load FMD ESD information
516 if (IsLoaded(kESD)) {
517 // AliInfo("Getting FMD event summary data");
518 Int_t read = fChainE->GetEntry(event);
519 if (read <= 0) return kFALSE;
520 fESD = fESDEvent->GetFMDData();
521 if (!fESD) return kFALSE;
522 }
523
524 // Possibly load FMD Digit information
525 if (IsLoaded(kRaw) || IsLoaded(kRawCalib)) {
526 Bool_t mon = fRawFile.Contains("mem://");
527 // AliInfo("Getting FMD raw data digits");
528 if (mon) std::cout << "Waiting for event ..." << std::flush;
529 do {
530 if (!fReader->NextEvent()) {
531 if (mon) {
532 gSystem->Sleep(3);
533 continue;
534 }
535 return kFALSE;
536 }
537 UInt_t eventType = fReader->GetType();
538 if(eventType == AliRawEventHeaderBase::kPhysicsEvent ||
539 eventType == AliRawEventHeaderBase::kCalibrationEvent)
540 break;
541 } while (true);
542 if (mon) std::cout << "got it" << std::endl;
543 // AliFMDRawReader r(fReader, 0);
544 fArrayA->Clear();
545 fFMDReader->ReadAdcs(fArrayA);
546 AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
547 }
548 fEventCount++;
549 return kTRUE;
550}
551
552
553//____________________________________________________________________
554Bool_t
555AliFMDInput::Event()
556{
557 // Process one event. The default implementation one or more of
558 //
559 // - ProcessHits if the hits are loaded.
560 // - ProcessDigits if the digits are loaded.
561 // - ProcessSDigits if the sumbable digits are loaded.
562 // - ProcessRecPoints if the reconstructed points are loaded.
563 // - ProcessESD if the event summary data is loaded
564 //
565 if (IsLoaded(kHits)) if (!ProcessHits()) return kFALSE;
566 if (IsLoaded(kTrackRefs))if (!ProcessTrackRefs()) return kFALSE;
567 if (IsLoaded(kKinematics) &&
568 IsLoaded(kHits)) if (!ProcessTracks()) return kFALSE;
569 if (IsLoaded(kKinematics))if (!ProcessStack()) return kFALSE;
570 if (IsLoaded(kSDigits)) if (!ProcessSDigits()) return kFALSE;
571 if (IsLoaded(kDigits)) if (!ProcessDigits()) return kFALSE;
572 if (IsLoaded(kRaw)) if (!ProcessRawDigits()) return kFALSE;
573 if (IsLoaded(kRawCalib)) if (!ProcessRawCalibDigits())return kFALSE;
574 if (IsLoaded(kRecPoints))if (!ProcessRecPoints()) return kFALSE;
575 if (IsLoaded(kESD)) if (!ProcessESDs()) return kFALSE;
576 if (IsLoaded(kUser)) if (!ProcessUsers()) return kFALSE;
577
578 return kTRUE;
579}
580
581//____________________________________________________________________
582Bool_t
583AliFMDInput::ProcessHits()
584{
585 // Read the hit tree, and pass each hit to the member function
586 // ProcessHit.
587 if (!fTreeH) {
588 AliError("No hit tree defined");
589 return kFALSE;
590 }
591 if (!fArrayH) {
592 AliError("No hit array defined");
593 return kFALSE;
594 }
595
596 Int_t nTracks = fTreeH->GetEntries();
597 for (Int_t i = 0; i < nTracks; i++) {
598 Int_t hitRead = fTreeH->GetEntry(i);
599 if (hitRead <= 0) continue;
600
601 Int_t nHit = fArrayH->GetEntries();
602 if (nHit <= 0) continue;
603
604 for (Int_t j = 0; j < nHit; j++) {
605 AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
606 if (!hit) continue;
607
608 TParticle* track = 0;
609 if (IsLoaded(kKinematics) && fStack) {
610 Int_t trackno = hit->Track();
611 track = fStack->Particle(trackno);
612 }
613 if (!ProcessHit(hit, track)) return kFALSE;
614 }
615 }
616 return kTRUE;
617}
618//____________________________________________________________________
619Bool_t
620AliFMDInput::ProcessTrackRefs()
621{
622 // Read the reconstrcted points tree, and pass each reconstruction
623 // object (AliFMDRecPoint) to either ProcessRecPoint.
624 if (!fTreeTR) {
625 AliError("No track reference tree defined");
626 return kFALSE;
627 }
628 if (!fArrayTR) {
629 AliError("No track reference array defined");
630 return kFALSE;
631 }
632
633 Int_t nEv = fTreeTR->GetEntries();
634 for (Int_t i = 0; i < nEv; i++) {
635 Int_t trRead = fTreeTR->GetEntry(i);
636 if (trRead <= 0) continue;
637 Int_t nTrackRefs = fArrayTR->GetEntries();
638 for (Int_t j = 0; j < nTrackRefs; j++) {
639 AliTrackReference* trackRef =
640 static_cast<AliTrackReference*>(fArrayTR->At(j));
641 if (!trackRef) continue;
642 // if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
643 TParticle* track = 0;
644 if (IsLoaded(kKinematics) && fStack) {
645 Int_t trackno = trackRef->GetTrack();
646 track = fStack->Particle(trackno);
647 }
648 if (!ProcessTrackRef(trackRef,track)) return kFALSE;
649 }
650 }
651 return kTRUE;
652}
653//____________________________________________________________________
654Bool_t
655AliFMDInput::ProcessTracks()
656{
657 // Read the hit tree, and pass each hit to the member function
658 // ProcessTrack.
659 if (!fStack) {
660 AliError("No track tree defined");
661 return kFALSE;
662 }
663 if (!fTreeH) {
664 AliError("No hit tree defined");
665 return kFALSE;
666 }
667 if (!fArrayH) {
668 AliError("No hit array defined");
669 return kFALSE;
670 }
671
672 // Int_t nTracks = fStack->GetNtrack();
673 Int_t nTracks = fTreeH->GetEntries();
674 for (Int_t i = 0; i < nTracks; i++) {
675 Int_t trackno = nTracks - i - 1;
676 TParticle* track = fStack->Particle(trackno);
677 if (!track) continue;
678
679 // Get the hits for this track.
680 Int_t hitRead = fTreeH->GetEntry(i);
681 Int_t nHit = fArrayH->GetEntries();
682 if (nHit == 0 || hitRead <= 0) {
683 // Let user code see the track, even if there's no hits.
684 if (!ProcessTrack(trackno, track, 0)) return kFALSE;
685 continue;
686 }
687
688 // Loop over the hits corresponding to this track.
689 for (Int_t j = 0; j < nHit; j++) {
690 AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
691 if (!ProcessTrack(trackno, track, hit)) return kFALSE;
692 }
693 }
694 return kTRUE;
695}
696//____________________________________________________________________
697Bool_t
698AliFMDInput::ProcessStack()
699{
700 // Read the hit tree, and pass each hit to the member function
701 // ProcessTrack.
702 if (!fStack) {
703 AliError("No track tree defined");
704 return kFALSE;
705 }
706 Int_t nTracks = fStack->GetNtrack();
707 for (Int_t i = 0; i < nTracks; i++) {
708 Int_t trackno = nTracks - i - 1;
709 TParticle* track = fStack->Particle(trackno);
710 if (!track) continue;
711
712 if (!ProcessParticle(trackno, track)) return kFALSE;
713 }
714 return kTRUE;
715}
716//____________________________________________________________________
717Bool_t
718AliFMDInput::ProcessDigits()
719{
720 // Read the digit tree, and pass each digit to the member function
721 // ProcessDigit.
722 if (!fTreeD) {
723 AliError("No digit tree defined");
724 return kFALSE;
725 }
726 if (!fArrayD) {
727 AliError("No digit array defined");
728 return kFALSE;
729 }
730
731 Int_t nEv = fTreeD->GetEntries();
732 for (Int_t i = 0; i < nEv; i++) {
733 Int_t digitRead = fTreeD->GetEntry(i);
734 if (digitRead <= 0) continue;
735 Int_t nDigit = fArrayD->GetEntries();
736 AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
737 if (nDigit <= 0) continue;
738 for (Int_t j = 0; j < nDigit; j++) {
739 AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
740 if (!digit) continue;
741 if (!ProcessDigit(digit)) return kFALSE;
742 }
743 }
744 return kTRUE;
745}
746
747//____________________________________________________________________
748Bool_t
749AliFMDInput::ProcessSDigits()
750{
751 // Read the summable digit tree, and pass each sumable digit to the
752 // member function ProcessSdigit.
753 if (!fTreeS) {
754 AliWarning("No sdigit tree defined");
755 return kTRUE; // Empty SDigits is fine
756 }
757 if (!fArrayS) {
758 AliWarning("No sdigit array defined");
759 return kTRUE; // Empty SDigits is fine
760 }
761
762 Int_t nEv = fTreeS->GetEntries();
763 for (Int_t i = 0; i < nEv; i++) {
764 Int_t sdigitRead = fTreeS->GetEntry(i);
765 if (sdigitRead <= 0) {
766 AliInfo(Form("Read nothing from tree"));
767 continue;
768 }
769 Int_t nSdigit = fArrayS->GetEntriesFast();
770 AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
771 AliInfo(Form("Got %5d digits for this event", nSdigit));
772 if (nSdigit <= 0) continue;
773 for (Int_t j = 0; j < nSdigit; j++) {
774 AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
775 if (!sdigit) continue;
776 if (!ProcessSDigit(sdigit)) return kFALSE;
777 }
778 }
779 return kTRUE;
780}
781
782//____________________________________________________________________
783Bool_t
784AliFMDInput::ProcessRawDigits()
785{
786 // Read the digit tree, and pass each digit to the member function
787 // ProcessDigit.
788 if (!fArrayA) {
789 AliError("No raw digit array defined");
790 return kFALSE;
791 }
792
793 Int_t nDigit = fArrayA->GetEntries();
794 if (nDigit <= 0) return kTRUE;
795 for (Int_t j = 0; j < nDigit; j++) {
796 AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
797 if (!digit) continue;
798 if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
799 digit->Print();
800 if (!ProcessRawDigit(digit)) return kFALSE;
801 }
802 return kTRUE;
803}
804
805//____________________________________________________________________
806Bool_t
807AliFMDInput::ProcessRawCalibDigits()
808{
809 // Read the digit tree, and pass each digit to the member function
810 // ProcessDigit.
811 if (!fArrayA) {
812 AliError("No raw digit array defined");
813 return kFALSE;
814 }
815
816 Int_t nDigit = fArrayA->GetEntries();
817 if (nDigit <= 0) return kTRUE;
818 for (Int_t j = 0; j < nDigit; j++) {
819 AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
820 if (!digit) continue;
821 if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
822 digit->Print();
823 if (!ProcessRawCalibDigit(digit)) return kFALSE;
824 }
825 return kTRUE;
826}
827
828//____________________________________________________________________
829Bool_t
830AliFMDInput::ProcessRecPoints()
831{
832 // Read the reconstrcted points tree, and pass each reconstruction
833 // object (AliFMDRecPoint) to either ProcessRecPoint.
834 if (!fTreeR) {
835 AliError("No recpoint tree defined");
836 return kFALSE;
837 }
838 if (!fArrayR) {
839 AliError("No recpoints array defined");
840 return kFALSE;
841 }
842
843 Int_t nEv = fTreeR->GetEntries();
844 for (Int_t i = 0; i < nEv; i++) {
845 Int_t recRead = fTreeR->GetEntry(i);
846 if (recRead <= 0) continue;
847 Int_t nRecPoint = fArrayR->GetEntries();
848 for (Int_t j = 0; j < nRecPoint; j++) {
849 AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
850 if (!recPoint) continue;
851 if (!ProcessRecPoint(recPoint)) return kFALSE;
852 }
853 }
854 return kTRUE;
855}
856
857//____________________________________________________________________
858Bool_t
859AliFMDInput::ProcessESDs()
860{
861 // Process event summary data
862 if (!fESD) return kFALSE;
863 for (UShort_t det = 1; det <= 3; det++) {
864 Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
865 for (Char_t* rng = rings; *rng != '\0'; rng++) {
866 UShort_t nsec = (*rng == 'I' ? 20 : 40);
867 UShort_t nstr = (*rng == 'I' ? 512 : 256);
868 for (UShort_t sec = 0; sec < nsec; sec++) {
869 for (UShort_t str = 0; str < nstr; str++) {
870 Float_t eta = fESD->Eta(det,*rng,sec,str);
871 Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
872 if (!fESD->IsAngleCorrected())
873 mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
874 if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
875 }
876 }
877 }
878 }
879 return kTRUE;
880}
881
882//____________________________________________________________________
883Bool_t
884AliFMDInput::ProcessUsers()
885{
886 // Process event summary data
887 for (UShort_t det = 1; det <= 3; det++) {
888 Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
889 for (Char_t* rng = rings; *rng != '\0'; rng++) {
890 UShort_t nsec = (*rng == 'I' ? 20 : 40);
891 UShort_t nstr = (*rng == 'I' ? 512 : 256);
892 for (UShort_t sec = 0; sec < nsec; sec++) {
893 for (UShort_t str = 0; str < nstr; str++) {
894 Float_t v = GetSignal(det,*rng,sec,str);
895 if (!ProcessUser(det, *rng, sec, str, v)) continue;
896 }
897 }
898 }
899 }
900 return kTRUE;
901}
902
903//____________________________________________________________________
904Bool_t
905AliFMDInput::End()
906{
907 // Called at the end of each event. Per default, it unloads the
908 // data trees and resets the pointers to the output arrays. Users
909 // can overload this, but should call this member function in the
910 // overloaded member function of the derived class.
911
912 // Check if we have been initialized
913 if (!fIsInit) {
914 AliError("Not initialized");
915 return fIsInit;
916 }
917 // Possibly unload global kinematics information
918 if (IsLoaded(kKinematics)) {
919 fLoader->UnloadKinematics();
920 // fTreeK = 0;
921 fStack = 0;
922 }
923 // Possibly unload FMD Hit information
924 if (IsLoaded(kHits)) {
925 fFMDLoader->UnloadHits();
926 fTreeH = 0;
927 }
928 // Possibly unload FMD Digit information
929 if (IsLoaded(kDigits)) {
930 fFMDLoader->UnloadDigits();
931 fTreeD = 0;
932 }
933 // Possibly unload FMD Sdigit information
934 if (IsLoaded(kSDigits)) {
935 fFMDLoader->UnloadSDigits();
936 fTreeS = 0;
937 }
938 // Possibly unload FMD RecPoints information
939 if (IsLoaded(kRecPoints)) {
940 fFMDLoader->UnloadRecPoints();
941 fTreeR = 0;
942 }
943 // AliInfo("Now out event");
944 return kTRUE;
945}
946
947//____________________________________________________________________
948Bool_t
949AliFMDInput::Run(UInt_t maxEvents)
950{
951 // Run over all events and files references in galice.root
952
953 Bool_t retval;
954 if (!(retval = Init())) return retval;
955
956 fNEvents = NEvents();
957 if (fNEvents < 0) fNEvents = maxEvents;
958 else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
959
960 Int_t event = 0;
961 for (; fNEvents < 0 || event < fNEvents; event++) {
962 printf("\rEvent %8d/%8d ...", event, fNEvents);
963 if (!(retval = Begin(event))) break;
964 if (!(retval = Event())) break;
965 if (!(retval = End())) break;
966 }
967 printf("Looped over %8d events\n", event+1);
968 if (!retval) return retval;
969 retval = Finish();
970 return retval;
971}
972
973//__________________________________________________________________
974TArrayF
975AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
976{
977 // Service function to define a logarithmic axis.
978 // Parameters:
979 // n Number of bins
980 // min Minimum of axis
981 // max Maximum of axis
982 TArrayF bins(n+1);
983 bins[0] = min;
984 if (n <= 20) {
985 for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
986 return bins;
987 }
988 Float_t dp = n / TMath::Log10(max / min);
989 Float_t pmin = TMath::Log10(min);
990 for (Int_t i = 1; i < n+1; i++) {
991 Float_t p = pmin + i / dp;
992 bins[i] = TMath::Power(10, p);
993 }
994 return bins;
995}
996
997//____________________________________________________________________
998void
999AliFMDInput::ScanDirectory(TSystemDirectory* dir,
1000 const TString& olddir,
1001 TChain* chain,
1002 const char* pattern, bool recursive)
1003{
1004 // Get list of files, and go back to old working directory
1005 TString oldDir(gSystem->WorkingDirectory());
1006 TList* files = dir->GetListOfFiles();
1007 gSystem->ChangeDirectory(oldDir);
1008
1009 // Sort list of files and check if we should add it
1010 if (!files) return;
1011 files->Sort();
1012 TIter next(files);
1013 TSystemFile* file = 0;
1014 while ((file = static_cast<TSystemFile*>(next()))) {
1015 TString name(file->GetName());
1016
1017 // Ignore special links
1018 if (name == "." || name == "..") continue;
1019
1020 // Check if this is a directory
1021 if (file->IsDirectory()) {
1022 if (recursive)
1023 ScanDirectory(static_cast<TSystemDirectory*>(file),
1024 olddir, chain,
1025 pattern,recursive);
1026 continue;
1027 }
1028
1029 // If this is not a root file, ignore
1030 if (!name.EndsWith(".root")) continue;
1031
1032 // If this file does not contain the pattern, ignore
1033 if (!name.Contains(pattern)) continue;
1034 if (name.Contains("friends")) continue;
1035
1036 // Get the path
1037 TString data(Form("%s/%s", file->GetTitle(), name.Data()));
1038
1039 TFile* test = TFile::Open(data.Data(), "READ");
1040 if (!test || test->IsZombie()) {
1041 ::Warning("ScanDirectory", "Failed to open file %s", data.Data());
1042 continue;
1043 }
1044 test->Close();
1045 chain->Add(data);
1046 }
1047}
1048
1049//____________________________________________________________________
1050TChain*
1051AliFMDInput::MakeChain(const char* what, const char* datadir, bool recursive)
1052{
1053 TString w(what);
1054 w.ToUpper();
1055 const char* treeName = 0;
1056 const char* pattern = 0;
1057 if (w.Contains("ESD")) { treeName = "esdTree"; pattern = "AliESD"; }
1058 else if (w.Contains("MC")) { treeName = "TE"; pattern = "galice"; }
1059 else {
1060 ::Error("MakeChain", "Unknown mode '%s' (not one of ESD, or MC)", what);
1061 return 0;
1062 }
1063
1064 // --- Our data chain ----------------------------------------------
1065 TChain* chain = new TChain(treeName);
1066
1067 // --- Get list of ESDs --------------------------------------------
1068 // Open source directory, and make sure we go back to were we were
1069 TString oldDir(gSystem->WorkingDirectory());
1070 TSystemDirectory d(datadir, datadir);
1071 ScanDirectory(&d, oldDir, chain, pattern, recursive);
1072
1073 return chain;
1074}
1075
1076
1077//____________________________________________________________________
1078//
1079// EOF
1080//