New RAW I/O. I rolled my own, because I wasn't happy with the old
[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
16 /* $Id$ */
17
18 //___________________________________________________________________
19 //
20 // The classes defined here, are utility classes for reading in data
21 // for the FMD.  They are  put in a seperate library to not polute the
22 // normal libraries.  The classes are intended to be used as base
23 // classes for customized class that do some sort of analysis on the
24 // various types of data produced by the FMD. 
25 //
26 // Latest changes by Christian Holm Christensen
27 //
28 #include "AliFMDInput.h"        // ALIFMDHIT_H
29 #include "AliLog.h"             // ALILOG_H
30 #include "AliLoader.h"          // ALILOADER_H
31 #include "AliRunLoader.h"       // ALIRUNLOADER_H
32 #include "AliRun.h"             // ALIRUN_H
33 #include "AliStack.h"           // ALISTACK_H
34 #include "AliFMD.h"             // ALIFMD_H
35 #include "AliFMDHit.h"          // ALIFMDHIT_H
36 #include "AliFMDDigit.h"        // ALIFMDDigit_H
37 #include "AliFMDRecPoint.h"     // ALIFMDRECPOINT_H
38 #include <AliESD.h>
39 #include <AliESDFMD.h>
40 #include <AliCDBManager.h>
41 #include <AliCDBEntry.h>
42 #include <AliAlignObjAngles.h>
43 #include <TTree.h>              // ROOT_TTree
44 #include <TChain.h>             // ROOT_TChain
45 #include <TParticle.h>          // ROOT_TParticle
46 #include <TString.h>            // ROOT_TString
47 #include <TDatabasePDG.h>       // ROOT_TDatabasePDG
48 #include <TMath.h>              // ROOT_TMath
49 #include <TGeoManager.h>        // ROOT_TGeoManager 
50 #include <TSystemDirectory.h>   // ROOT_TSystemDirectory
51 #include <Riostream.h>          // ROOT_Riostream
52
53 //____________________________________________________________________
54 ClassImp(AliFMDInput)
55 #if 0
56   ; // This is here to keep Emacs for indenting the next line
57 #endif
58
59
60 //____________________________________________________________________
61 AliFMDInput::AliFMDInput()
62   : fGAliceFile(""), 
63     fLoader(0),
64     fRun(0), 
65     fStack(0),
66     fFMDLoader(0), 
67     fFMD(0),
68     fMainESD(0),
69     fESD(0),
70     fTreeE(0),
71     fTreeH(0),
72     fTreeD(0),
73     fTreeS(0),
74     fTreeR(0), 
75     fChainE(0),
76     fArrayE(0),
77     fArrayH(0),
78     fArrayD(0),
79     fArrayS(0), 
80     fArrayR(0), 
81     fTreeMask(0), 
82     fIsInit(kFALSE)
83 {
84   // Constructor of an FMD input object.  Specify what data to read in
85   // using the AddLoad member function.   Sub-classes should at a
86   // minimum overload the member function Event.   A full job can be
87   // executed using the member function Run. 
88 }
89
90   
91
92 //____________________________________________________________________
93 AliFMDInput::AliFMDInput(const char* gAliceFile)
94   : fGAliceFile(gAliceFile),
95     fLoader(0),
96     fRun(0), 
97     fStack(0),
98     fFMDLoader(0), 
99     fFMD(0),
100     fMainESD(0),
101     fESD(0),
102     fTreeE(0),
103     fTreeH(0),
104     fTreeD(0),
105     fTreeS(0),
106     fTreeR(0), 
107     fChainE(0),
108     fArrayE(0),
109     fArrayH(0),
110     fArrayD(0),
111     fArrayS(0), 
112     fArrayR(0), 
113     fTreeMask(0), 
114     fIsInit(kFALSE)
115 {
116   // Constructor of an FMD input object.  Specify what data to read in
117   // using the AddLoad member function.   Sub-classes should at a
118   // minimum overload the member function Event.   A full job can be
119   // executed using the member function Run. 
120 }
121
122 //____________________________________________________________________
123 Int_t
124 AliFMDInput::NEvents() const 
125 {
126   // Get number of events
127   if (fTreeE) return fTreeE->GetEntries();
128   return -1;
129 }
130
131 //____________________________________________________________________
132 Bool_t
133 AliFMDInput::Init()
134 {
135   // Initialize the object.  Get the needed loaders, and such. 
136
137   // Check if we have been initialized
138   if (fIsInit) { 
139     AliWarning("Already initialized");
140     return fIsInit;
141   }
142   if (fGAliceFile.IsNull()) fGAliceFile = "galice.root";
143   // Get the loader
144   fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
145   if (!fLoader) {
146     AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
147     return kFALSE;
148   }
149   
150   // Get the run 
151   if  (fLoader->LoadgAlice()) return kFALSE;
152   fRun = fLoader->GetAliRun();
153   
154   // Get the FMD 
155   fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
156   if (!fFMD) {
157     AliError("Failed to get detector FMD from loader");
158     return kFALSE;
159   }
160   
161   // Get the FMD loader
162   fFMDLoader = fLoader->GetLoader("FMDLoader");
163   if (!fFMDLoader) {
164     AliError("Failed to get detector FMD loader from loader");
165     return kFALSE;
166   }
167   if (fLoader->LoadHeader()) { 
168     AliError("Failed to get event header information from loader");
169     return kFALSE;
170   }
171   fTreeE = fLoader->TreeE();
172
173   // Optionally, get the ESD files
174   if (TESTBIT(fTreeMask, kESD)) {
175     fChainE = new TChain("esdTree");
176     TSystemDirectory dir;
177     TList*           files = dir.GetListOfFiles();
178     TSystemFile*     file;
179     files->Sort();
180     TIter            next(files);
181     while ((file = static_cast<TSystemFile*>(next()))) {
182       TString fname(file->GetName());
183       if (fname.Contains("AliESDs")) fChainE->AddFile(fname.Data());
184     }
185     fChainE->SetBranchAddress("ESD", &fMainESD);
186   }
187     
188   // Optionally, get the geometry 
189   if (TESTBIT(fTreeMask, kGeometry)) {
190     TString fname(fRun->GetGeometryFileName());
191     if (fname.IsNull()) {
192       Warning("Init", "No file name for the geometry from AliRun");
193       fname = gSystem->DirName(fGAliceFile);
194       fname.Append("/geometry.root");
195     }
196     fGeoManager = TGeoManager::Import(fname.Data());
197     if (!fGeoManager) {
198       Fatal("Init", "No geometry manager found");
199       return kFALSE;
200     }
201     AliCDBManager* cdb   = AliCDBManager::Instance();
202     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
203     if (align) {
204       AliInfo("Got alignment data from CDB");
205       TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
206       if (!array) {
207         AliWarning("Invalid align data from CDB");
208       }
209       else {
210         Int_t nAlign = array->GetEntries();
211         for (Int_t i = 0; i < nAlign; i++) {
212           AliAlignObjAngles* a = static_cast<AliAlignObjAngles*>(array->At(i));
213           if (!a->ApplyToGeometry()) {
214             AliWarning(Form("Failed to apply alignment to %s", 
215                             a->GetVolPath()));
216           }
217         }
218       }
219     }
220   }
221
222   
223   fIsInit = kTRUE;
224   return fIsInit;
225 }
226
227 //____________________________________________________________________
228 Bool_t
229 AliFMDInput::Begin(Int_t event)
230 {
231   // Called at the begining of each event.  Per default, it gets the
232   // data trees and gets pointers to the output arrays.   Users can
233   // overload this, but should call this member function in the
234   // overloaded member function of the derived class. 
235
236   // Check if we have been initialized
237   if (!fIsInit) { 
238     AliError("Not initialized");
239     return fIsInit;
240   }
241   // Get the event 
242   if (fLoader->GetEvent(event)) return kFALSE;
243   AliInfo(Form("Now in event %d/%d", event, NEvents()));
244
245   // Possibly load global kinematics information 
246   if (TESTBIT(fTreeMask, kKinematics)) {
247     AliInfo("Getting kinematics");
248     if (fLoader->LoadKinematics()) return kFALSE;
249     fStack = fLoader->Stack();
250   }
251   // Possibly load FMD Hit information 
252   if (TESTBIT(fTreeMask, kHits)) {
253     AliInfo("Getting FMD hits");
254     if (fFMDLoader->LoadHits()) return kFALSE;
255     fTreeH = fFMDLoader->TreeH();
256     if (!fArrayH) fArrayH = fFMD->Hits(); 
257   }
258   // Possibly load FMD Digit information 
259   if (TESTBIT(fTreeMask, kDigits)) {
260     AliInfo("Getting FMD digits");
261     if (fFMDLoader->LoadDigits()) return kFALSE;
262     fTreeD = fFMDLoader->TreeD();
263     if (!fArrayD) fArrayD = fFMD->Digits();
264   }
265   // Possibly load FMD Sdigit information 
266   if (TESTBIT(fTreeMask, kSDigits)) {
267     AliInfo("Getting FMD summable digits");
268     if (fFMDLoader->LoadSDigits()) return kFALSE;
269     fTreeS = fFMDLoader->TreeS();
270     if (!fArrayS) fArrayS = fFMD->SDigits();
271   }
272   // Possibly load FMD RecPoints information 
273   if (TESTBIT(fTreeMask, kRecPoints)) {
274     AliInfo("Getting FMD reconstructed points");
275     if (fFMDLoader->LoadRecPoints()) return kFALSE;
276     fTreeR = fFMDLoader->TreeR();
277     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
278     fTreeR->SetBranchAddress("FMD",  &fArrayR);
279   }
280
281   // Possibly load FMD ESD information 
282   if (TESTBIT(fTreeMask, kESD)) {
283     AliInfo("Getting FMD event summary data");
284     Int_t read = fChainE->GetEntry(event);
285     if (read <= 0) return kFALSE;
286     fESD = fMainESD->GetFMDData();
287     if (!fESD) return kFALSE;
288   }
289   
290   return kTRUE;
291 }
292
293
294 //____________________________________________________________________
295 Bool_t 
296 AliFMDInput::Event()
297 {
298   // Process one event.  The default implementation one or more of 
299   //
300   //   -  ProcessHits       if the hits are loaded. 
301   //   -  ProcessDigits     if the digits are loaded. 
302   //   -  ProcessSDigits    if the sumbable digits are loaded. 
303   //   -  ProcessRecPoints  if the reconstructed points are loaded. 
304   //   -  ProcessESD        if the event summary data is loaded
305   // 
306   if (TESTBIT(fTreeMask, kHits)) 
307     if (!ProcessHits()) return kFALSE; 
308   if (TESTBIT(fTreeMask, kDigits)) 
309     if (!ProcessDigits()) return kFALSE;
310   if (TESTBIT(fTreeMask, kSDigits)) 
311     if (!ProcessSDigits()) return kFALSE;
312   if (TESTBIT(fTreeMask, kRecPoints)) 
313     if (!ProcessRecPoints()) return kFALSE;
314   if (TESTBIT(fTreeMask, kESD))
315     if (!ProcessESD(fESD)) return kFALSE;
316   
317   return kTRUE;
318 }
319
320 //____________________________________________________________________
321 Bool_t 
322 AliFMDInput::ProcessHits()
323 {
324   if (!fTreeH) {
325     AliError("No hit tree defined");
326     return kFALSE;
327   }
328   Int_t nTracks = fTreeH->GetEntries();
329   for (Int_t i = 0; i < nTracks; i++) {
330     Int_t hitRead  = fTreeH->GetEntry(i);
331     if (hitRead <= 0) continue;
332     if (!fArrayH) {
333       AliError("No hit array defined");
334       return kFALSE;
335     }
336     Int_t nHit = fArrayH->GetEntries();
337     if (nHit <= 0) continue;
338     for (Int_t j = 0; j < nHit; j++) {
339       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
340       if (!hit) continue;
341       TParticle* track = 0;
342       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
343         Int_t trackno = hit->Track();
344         track = fStack->Particle(trackno);
345       }
346       if (!ProcessHit(hit, track)) return kFALSE;
347     }    
348   }
349   return kTRUE;
350 }
351
352 //____________________________________________________________________
353 Bool_t 
354 AliFMDInput::ProcessDigits()
355 {
356   // Read the digit tree, and pass each digit to the member function
357   // ProcessDigit.
358   Int_t nEv = fTreeD->GetEntries();
359   for (Int_t i = 0; i < nEv; i++) {
360     Int_t digitRead  = fTreeD->GetEntry(i);
361     if (digitRead <= 0) continue;
362     Int_t nDigit = fArrayD->GetEntries();
363     if (nDigit <= 0) continue;
364     for (Int_t j = 0; j < nDigit; j++) {
365       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
366       if (!digit) continue;
367       if (!ProcessDigit(digit)) return kFALSE;
368     }    
369   }
370   return kTRUE;
371 }
372
373 //____________________________________________________________________
374 Bool_t 
375 AliFMDInput::ProcessSDigits()
376 {
377   // Read the summable digit tree, and pass each sumable digit to the
378   // member function ProcessSdigit.
379   Int_t nEv = fTreeD->GetEntries();
380   for (Int_t i = 0; i < nEv; i++) {
381     Int_t sdigitRead  = fTreeS->GetEntry(i);
382     if (sdigitRead <= 0) continue;
383     Int_t nSdigit = fArrayS->GetEntries();
384     if (nSdigit <= 0) continue;
385     for (Int_t j = 0; j < nSdigit; j++) {
386       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
387       if (!sdigit) continue;
388       if (!ProcessSDigit(sdigit)) return kFALSE;
389     }    
390   }
391   return kTRUE;
392 }
393
394 //____________________________________________________________________
395 Bool_t 
396 AliFMDInput::ProcessRecPoints()
397 {
398   // Read the reconstrcted points tree, and pass each reconstruction
399   // object (AliFMDRecPoint) to either ProcessRecPoint.
400   Int_t nEv = fTreeR->GetEntries();
401   for (Int_t i = 0; i < nEv; i++) {
402     Int_t recRead  = fTreeR->GetEntry(i);
403     if (recRead <= 0) continue;
404     Int_t nRecPoint = fArrayR->GetEntries();
405     for (Int_t j = 0; j < nRecPoint; j++) {
406       AliFMDRecPoint* recPoint = static_cast<AliFMDRecPoint*>(fArrayR->At(j));
407       if (!recPoint) continue;
408       if (!ProcessRecPoint(recPoint)) return kFALSE;
409     }    
410   }
411   return kTRUE;
412 }
413
414 //____________________________________________________________________
415 Bool_t
416 AliFMDInput::End()
417 {
418   // Called at the end of each event.  Per default, it unloads the
419   // data trees and resets the pointers to the output arrays.   Users
420   // can overload this, but should call this member function in the
421   // overloaded member function of the derived class. 
422
423   // Check if we have been initialized
424   if (!fIsInit) { 
425     AliError("Not initialized");
426     return fIsInit;
427   }
428   // Possibly unload global kinematics information 
429   if (TESTBIT(fTreeMask, kKinematics)) {
430     fLoader->UnloadKinematics();
431     // fTreeK = 0;
432     fStack = 0;
433   }
434   // Possibly unload FMD Hit information 
435   if (TESTBIT(fTreeMask, kHits)) {
436     fFMDLoader->UnloadHits();
437     fTreeH = 0;
438   }
439   // Possibly unload FMD Digit information 
440   if (TESTBIT(fTreeMask, kDigits)) {
441     fFMDLoader->UnloadDigits();
442     fTreeD = 0;
443   }
444   // Possibly unload FMD Sdigit information 
445   if (TESTBIT(fTreeMask, kSDigits)) {
446     fFMDLoader->UnloadSDigits();
447     fTreeS = 0;
448   }
449   // Possibly unload FMD RecPoints information 
450   if (TESTBIT(fTreeMask, kRecPoints)) {
451     fFMDLoader->UnloadRecPoints();
452     fTreeR = 0;
453   }
454   AliInfo("Now out event");
455   return kTRUE;
456 }
457
458 //____________________________________________________________________
459 Bool_t
460 AliFMDInput::Run()
461 {
462   // Run over all events and files references in galice.root 
463
464   Bool_t retval;
465   if (!(retval = Init())) return retval;
466
467   Int_t nEvents = NEvents();
468   for (Int_t event = 0; event < nEvents; event++) {
469     if (!(retval = Begin(event))) break;
470     if (!(retval = Event())) break;
471     if (!(retval = End())) break;
472   }
473   if (!retval) return retval;
474   retval = Finish();
475   return retval;
476 }
477
478 //====================================================================
479 ClassImp(AliFMDInputHits)
480 #if 0
481   ;
482 #endif
483
484
485 //====================================================================
486 ClassImp(AliFMDInputDigits)
487 #if 0
488   ;
489 #endif
490
491 //====================================================================
492 ClassImp(AliFMDInputSDigits)
493 #if 0
494   ;
495 #endif
496
497 //====================================================================
498 ClassImp(AliFMDInputRecPoints)
499 #if 0
500   ;
501 #endif
502
503 //____________________________________________________________________
504 //
505 // EOF
506 //