1 /**************************************************************************
2 * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
21 //___________________________________________________________________
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.
29 // Latest changes by Christian Holm Christensen
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
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>
66 //____________________________________________________________________
69 ; // This is here to keep Emacs for indenting the next line
73 //____________________________________________________________________
74 AliFMDInput::AliFMDInput()
75 : TNamed("AliFMDInput", "Input handler for various FMD data"),
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.
117 //____________________________________________________________________
118 AliFMDInput::AliFMDInput(const char* gAliceFile)
119 : TNamed("AliFMDInput", "Input handler for various FMD data"),
120 fGAliceFile(gAliceFile),
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.
159 //____________________________________________________________________
161 AliFMDInput::NEvents() const
163 // Get number of events
164 if (TESTBIT(fTreeMask, kRaw)) return fReader->GetNumberOfEvents();
165 if (fChainE) return fChainE->GetEntriesFast();
166 if (fTreeE) return fTreeE->GetEntries();
170 //____________________________________________________________________
174 // Initialize the object. Get the needed loaders, and such.
176 // Check if we have been initialized
178 AliWarning("Already initialized");
181 Info("Init","Initialising w/mask 0x%04x\n"
194 TESTBIT(fTreeMask, kHits),
195 TESTBIT(fTreeMask, kKinematics),
196 TESTBIT(fTreeMask, kDigits),
197 TESTBIT(fTreeMask, kSDigits),
198 TESTBIT(fTreeMask, kHeader),
199 TESTBIT(fTreeMask, kRecPoints),
200 TESTBIT(fTreeMask, kESD),
201 TESTBIT(fTreeMask, kRaw),
202 TESTBIT(fTreeMask, kGeometry),
203 TESTBIT(fTreeMask, kTracks),
204 TESTBIT(fTreeMask, kTrackRefs));
206 if (TESTBIT(fTreeMask, kDigits) ||
207 TESTBIT(fTreeMask, kSDigits) ||
208 TESTBIT(fTreeMask, kKinematics) ||
209 TESTBIT(fTreeMask, kTrackRefs) ||
210 TESTBIT(fTreeMask, kTracks) ||
211 TESTBIT(fTreeMask, kHeader)) {
212 if (!gSystem->FindFile(".:/", fGAliceFile)) {
213 AliWarning(Form("Cannot find file %s in .:/", fGAliceFile.Data()));
216 fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
218 AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
221 AliInfo(Form("Opened GAlice file %s", fGAliceFile.Data()));
223 if (fLoader->LoadgAlice()) return kFALSE;
225 fRun = fLoader->GetAliRun();
228 fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
230 AliError("Failed to get detector FMD from loader");
234 // Get the FMD loader
235 fFMDLoader = fLoader->GetLoader("FMDLoader");
237 AliError("Failed to get detector FMD loader from loader");
240 if (fLoader->LoadHeader()) {
241 AliError("Failed to get event header information from loader");
244 fTreeE = fLoader->TreeE();
248 // Optionally, get the ESD files
249 if (TESTBIT(fTreeMask, kESD)) {
250 fChainE = new TChain("esdTree");
251 TSystemDirectory dir(".",".");
252 TList* files = dir.GetListOfFiles();
253 TSystemFile* file = 0;
255 AliError("No files");
260 while ((file = static_cast<TSystemFile*>(next()))) {
261 TString fname(file->GetName());
262 if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
264 fESDEvent = new AliESDEvent();
265 fESDEvent->ReadFromTree(fChainE);
266 // fChainE->SetBranchAddress("ESD", &fMainESD);
270 if (TESTBIT(fTreeMask, kRaw)) {
271 AliInfo("Getting FMD raw data digits");
272 fArrayA = new TClonesArray("AliFMDDigit");
274 if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
275 fReader = new AliRawReaderRoot(fRawFile.Data());
276 else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
277 fReader = new AliRawReaderDate(fRawFile.Data());
279 fReader = new AliRawReaderFile(-1);
281 if(!fRawFile.IsNull())
282 fReader = AliRawReader::Create(fRawFile.Data());
284 fReader = new AliRawReaderFile(-1);
286 fFMDReader = new AliFMDRawReader(fReader, 0);
289 // Optionally, get the geometry
290 if (TESTBIT(fTreeMask, kGeometry)) {
292 TString fname(fRun->GetGeometryFileName());
293 if (fname.IsNull()) {
294 Warning("Init", "No file name for the geometry from AliRun");
295 fname = gSystem->DirName(fGAliceFile);
296 fname.Append("/geometry.root");
298 fGeoManager = TGeoManager::Import(fname.Data());
300 Fatal("Init", "No geometry manager found");
305 AliGeomManager::LoadGeometry();
307 AliCDBManager* cdb = AliCDBManager::Instance();
308 AliCDBEntry* align = cdb->Get("FMD/Align/Data");
310 AliInfo("Got alignment data from CDB");
311 TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
313 AliWarning("Invalid align data from CDB");
316 Int_t nAlign = array->GetEntries();
317 for (Int_t i = 0; i < nAlign; i++) {
318 AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
319 if (!a->ApplyToGeometry()) {
320 AliWarning(Form("Failed to apply alignment to %s",
333 //____________________________________________________________________
335 AliFMDInput::Begin(Int_t event)
337 // Called at the begining of each event. Per default, it gets the
338 // data trees and gets pointers to the output arrays. Users can
339 // overload this, but should call this member function in the
340 // overloaded member function of the derived class.
342 // Check if we have been initialized
344 AliError("Not initialized");
349 if (fLoader && fLoader->GetEvent(event)) return kFALSE;
350 AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
352 // Possibly load global kinematics information
353 if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
354 // AliInfo("Getting kinematics");
355 if (fLoader->LoadKinematics("READ")) return kFALSE;
356 fStack = fLoader->Stack();
359 // Possibly load FMD Hit information
360 if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
361 // AliInfo("Getting FMD hits");
362 if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
363 fTreeH = fFMDLoader->TreeH();
364 if (!fArrayH) fArrayH = fFMD->Hits();
367 // Possibly load FMD TrackReference information
368 if (TESTBIT(fTreeMask, kTrackRefs) || TESTBIT(fTreeMask, kTracks)) {
369 // AliInfo("Getting FMD hits");
370 if (!fLoader || fLoader->LoadTrackRefs("READ")) return kFALSE;
371 fTreeTR = fLoader->TreeTR();
372 if (!fArrayTR) fArrayTR = new TClonesArray("AliTrackReference");
373 fTreeTR->SetBranchAddress("TrackReferences", &fArrayTR);
376 // Possibly load heaedr information
377 if (TESTBIT(fTreeMask, kHeader)) {
378 // AliInfo("Getting FMD hits");
379 if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
380 fHeader = fLoader->GetHeader();
383 // Possibly load FMD Digit information
384 if (TESTBIT(fTreeMask, kDigits)) {
385 // AliInfo("Getting FMD digits");
386 if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
387 fTreeD = fFMDLoader->TreeD();
389 if (!fArrayD) fArrayD = fFMD->Digits();
393 AliWarning(Form("Failed to load FMD Digits"));
397 // Possibly load FMD Sdigit information
398 if (TESTBIT(fTreeMask, kSDigits)) {
399 // AliInfo("Getting FMD summable digits");
400 if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) {
401 AliWarning("Failed to load SDigits!");
404 fTreeS = fFMDLoader->TreeS();
405 if (!fArrayS) fArrayS = fFMD->SDigits();
408 // Possibly load FMD RecPoints information
409 if (TESTBIT(fTreeMask, kRecPoints)) {
410 // AliInfo("Getting FMD reconstructed points");
411 if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
412 fTreeR = fFMDLoader->TreeR();
413 if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
414 fTreeR->SetBranchAddress("FMD", &fArrayR);
417 // Possibly load FMD ESD information
418 if (TESTBIT(fTreeMask, kESD)) {
419 // AliInfo("Getting FMD event summary data");
420 Int_t read = fChainE->GetEntry(event);
421 if (read <= 0) return kFALSE;
422 fESD = fESDEvent->GetFMDData();
423 if (!fESD) return kFALSE;
425 TFile* f = fChainE->GetFile();
427 TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
429 TStreamerInfo* info = static_cast<TStreamerInfo*>(o);
430 std::cout << "AliFMDMap class version read is "
431 << info->GetClassVersion() << std::endl;
434 // fESD->CheckNeedUShort(fChainE->GetFile());
438 // Possibly load FMD Digit information
439 if (TESTBIT(fTreeMask, kRaw)) {
440 // AliInfo("Getting FMD raw data digits");
441 if (!fReader->NextEvent()) return kFALSE;
442 // AliFMDRawReader r(fReader, 0);
444 fFMDReader->ReadAdcs(fArrayA);
445 AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
452 //____________________________________________________________________
456 // Process one event. The default implementation one or more of
458 // - ProcessHits if the hits are loaded.
459 // - ProcessDigits if the digits are loaded.
460 // - ProcessSDigits if the sumbable digits are loaded.
461 // - ProcessRecPoints if the reconstructed points are loaded.
462 // - ProcessESD if the event summary data is loaded
464 if (TESTBIT(fTreeMask, kHits))
465 if (!ProcessHits()) return kFALSE;
466 if (TESTBIT(fTreeMask, kTrackRefs))
467 if (!ProcessTrackRefs()) return kFALSE;
468 if (TESTBIT(fTreeMask, kTracks))
469 if (!ProcessTracks()) return kFALSE;
470 if (TESTBIT(fTreeMask, kSDigits))
471 if (!ProcessSDigits()) return kFALSE;
472 if (TESTBIT(fTreeMask, kDigits))
473 if (!ProcessDigits()) return kFALSE;
474 if (TESTBIT(fTreeMask, kRaw))
475 if (!ProcessRawDigits()) return kFALSE;
476 if (TESTBIT(fTreeMask, kRecPoints))
477 if (!ProcessRecPoints()) return kFALSE;
478 if (TESTBIT(fTreeMask, kESD))
479 if (!ProcessESDs()) return kFALSE;
484 //____________________________________________________________________
486 AliFMDInput::ProcessHits()
488 // Read the hit tree, and pass each hit to the member function
491 AliError("No hit tree defined");
495 AliError("No hit array defined");
499 Int_t nTracks = fTreeH->GetEntries();
500 for (Int_t i = 0; i < nTracks; i++) {
501 Int_t hitRead = fTreeH->GetEntry(i);
502 if (hitRead <= 0) continue;
504 Int_t nHit = fArrayH->GetEntries();
505 if (nHit <= 0) continue;
507 for (Int_t j = 0; j < nHit; j++) {
508 AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
511 TParticle* track = 0;
512 if (TESTBIT(fTreeMask, kKinematics) && fStack) {
513 Int_t trackno = hit->Track();
514 track = fStack->Particle(trackno);
516 if (!ProcessHit(hit, track)) return kFALSE;
521 //____________________________________________________________________
523 AliFMDInput::ProcessTrackRefs()
525 // Read the reconstrcted points tree, and pass each reconstruction
526 // object (AliFMDRecPoint) to either ProcessRecPoint.
528 AliError("No track reference tree defined");
532 AliError("No track reference array defined");
536 Int_t nEv = fTreeTR->GetEntries();
537 for (Int_t i = 0; i < nEv; i++) {
538 Int_t trRead = fTreeTR->GetEntry(i);
539 if (trRead <= 0) continue;
540 Int_t nTrackRefs = fArrayTR->GetEntries();
541 for (Int_t j = 0; j < nTrackRefs; j++) {
542 AliTrackReference* trackRef = static_cast<AliTrackReference*>(fArrayTR->At(j));
543 if (!trackRef) continue;
544 if (trackRef->DetectorId() != AliTrackReference::kFMD) continue;
545 TParticle* track = 0;
546 if (TESTBIT(fTreeMask, kKinematics) && fStack) {
547 Int_t trackno = trackRef->GetTrack();
548 track = fStack->Particle(trackno);
550 if (!ProcessTrackRef(trackRef,track)) return kFALSE;
555 //____________________________________________________________________
557 AliFMDInput::ProcessTracks()
559 // Read the hit tree, and pass each hit to the member function
562 AliError("No track tree defined");
566 AliError("No hit tree defined");
570 AliError("No hit array defined");
574 // Int_t nTracks = fStack->GetNtrack();
575 Int_t nTracks = fTreeH->GetEntries();
576 for (Int_t i = 0; i < nTracks; i++) {
577 Int_t trackno = nTracks - i - 1;
578 TParticle* track = fStack->Particle(trackno);
579 if (!track) continue;
581 // Get the hits for this track.
582 Int_t hitRead = fTreeH->GetEntry(i);
583 Int_t nHit = fArrayH->GetEntries();
584 if (nHit == 0 || hitRead <= 0) {
585 // Let user code see the track, even if there's no hits.
586 if (!ProcessTrack(trackno, track, 0)) return kFALSE;
590 // Loop over the hits corresponding to this track.
591 for (Int_t j = 0; j < nHit; j++) {
592 AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
593 if (!ProcessTrack(trackno, track, hit)) return kFALSE;
599 //____________________________________________________________________
601 AliFMDInput::ProcessDigits()
603 // Read the digit tree, and pass each digit to the member function
606 AliError("No digit tree defined");
610 AliError("No digit array defined");
614 Int_t nEv = fTreeD->GetEntries();
615 for (Int_t i = 0; i < nEv; i++) {
616 Int_t digitRead = fTreeD->GetEntry(i);
617 if (digitRead <= 0) continue;
618 Int_t nDigit = fArrayD->GetEntries();
619 AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
620 if (nDigit <= 0) continue;
621 for (Int_t j = 0; j < nDigit; j++) {
622 AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
623 if (!digit) continue;
624 if (!ProcessDigit(digit)) return kFALSE;
630 //____________________________________________________________________
632 AliFMDInput::ProcessSDigits()
634 // Read the summable digit tree, and pass each sumable digit to the
635 // member function ProcessSdigit.
637 AliWarning("No sdigit tree defined");
638 return kTRUE; // Empty SDigits is fine
641 AliWarning("No sdigit array defined");
642 return kTRUE; // Empty SDigits is fine
645 Int_t nEv = fTreeS->GetEntries();
646 for (Int_t i = 0; i < nEv; i++) {
647 Int_t sdigitRead = fTreeS->GetEntry(i);
648 if (sdigitRead <= 0) {
649 AliInfo(Form("Read nothing from tree"));
652 Int_t nSdigit = fArrayS->GetEntriesFast();
653 AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
654 AliInfo(Form("Got %5d digits for this event", nSdigit));
655 if (nSdigit <= 0) continue;
656 for (Int_t j = 0; j < nSdigit; j++) {
657 AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
658 if (!sdigit) continue;
659 if (!ProcessSDigit(sdigit)) return kFALSE;
665 //____________________________________________________________________
667 AliFMDInput::ProcessRawDigits()
669 // Read the digit tree, and pass each digit to the member function
672 AliError("No raw digit array defined");
676 Int_t nDigit = fArrayA->GetEntries();
677 if (nDigit <= 0) return kTRUE;
678 for (Int_t j = 0; j < nDigit; j++) {
679 AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
680 if (!digit) continue;
681 if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30)
683 if (!ProcessRawDigit(digit)) return kFALSE;
688 //____________________________________________________________________
690 AliFMDInput::ProcessRecPoints()
692 // Read the reconstrcted points tree, and pass each reconstruction
693 // object (AliFMDRecPoint) to either ProcessRecPoint.
695 AliError("No recpoint tree defined");
699 AliError("No recpoints array defined");
703 Int_t nEv = fTreeR->GetEntries();
704 for (Int_t i = 0; i < nEv; i++) {
705 Int_t recRead = fTreeR->GetEntry(i);
706 if (recRead <= 0) continue;
707 Int_t nRecPoint = fArrayR->GetEntries();
708 for (Int_t j = 0; j < nRecPoint; j++) {
709 AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
710 if (!recPoint) continue;
711 if (!ProcessRecPoint(recPoint)) return kFALSE;
717 //____________________________________________________________________
719 AliFMDInput::ProcessESDs()
721 // Process event summary data
722 if (!fESD) return kFALSE;
723 for (UShort_t det = 1; det <= 3; det++) {
724 Char_t rings[] = { 'I', (det == 1 ? '\0' : 'O'), '\0' };
725 for (Char_t* rng = rings; *rng != '\0'; rng++) {
726 UShort_t nsec = (*rng == 'I' ? 20 : 40);
727 UShort_t nstr = (*rng == 'I' ? 512 : 256);
728 for (UShort_t sec = 0; sec < nsec; sec++) {
729 for (UShort_t str = 0; str < nstr; str++) {
730 Float_t eta = fESD->Eta(det,*rng,sec,str);
731 Float_t mult = fESD->Multiplicity(det,*rng,sec,str);
732 if (!fESD->IsAngleCorrected())
733 mult *= TMath::Abs(TMath::Cos(2.*TMath::ATan(TMath::Exp(-eta))));
734 if (!ProcessESD(det, *rng, sec, str, eta, mult)) continue;
742 //____________________________________________________________________
746 // Called at the end of each event. Per default, it unloads the
747 // data trees and resets the pointers to the output arrays. Users
748 // can overload this, but should call this member function in the
749 // overloaded member function of the derived class.
751 // Check if we have been initialized
753 AliError("Not initialized");
756 // Possibly unload global kinematics information
757 if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
758 fLoader->UnloadKinematics();
762 // Possibly unload FMD Hit information
763 if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
764 fFMDLoader->UnloadHits();
767 // Possibly unload FMD Digit information
768 if (TESTBIT(fTreeMask, kDigits)) {
769 fFMDLoader->UnloadDigits();
772 // Possibly unload FMD Sdigit information
773 if (TESTBIT(fTreeMask, kSDigits)) {
774 fFMDLoader->UnloadSDigits();
777 // Possibly unload FMD RecPoints information
778 if (TESTBIT(fTreeMask, kRecPoints)) {
779 fFMDLoader->UnloadRecPoints();
782 // AliInfo("Now out event");
786 //____________________________________________________________________
790 // Run over all events and files references in galice.root
793 if (!(retval = Init())) return retval;
795 Int_t nEvents = NEvents();
796 for (Int_t event = 0; nEvents < 0 || event < nEvents; event++) {
797 if (!(retval = Begin(event))) break;
798 if (!(retval = Event())) break;
799 if (!(retval = End())) break;
801 if (!retval) return retval;
806 //__________________________________________________________________
808 AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max)
810 // Service function to define a logarithmic axis.
813 // min Minimum of axis
814 // max Maximum of axis
818 for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
821 Float_t dp = n / TMath::Log10(max / min);
822 Float_t pmin = TMath::Log10(min);
823 for (Int_t i = 1; i < n+1; i++) {
824 Float_t p = pmin + i / dp;
825 bins[i] = TMath::Power(10, p);
832 //____________________________________________________________________