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