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