13-oct-2005 NvE TTask derived class IceCleanHits introduced to perform hit cleaning.
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 // This is a class that constructs AliFMDMult (reconstructed
21 // multiplicity) from of Digits
22 //
23 // This class reads either digits from a TClonesArray or raw data from
24 // a DDL file (or similar), and stores the read ADC counts in an
25 // internal cache (fAdcs). 
26 //
27 // From the cached values it then calculates the number of particles
28 // that hit a region of the FMDs, as specified by the user. 
29 //
30 // The reconstruction can be done in two ways: Either via counting the
31 // number of empty strips (Poisson method), or by converting the ADC
32 // signal to an energy deposition, and then dividing by the typical
33 // energy loss of a particle.
34 // 
35 //      +---------------------+       +---------------------+
36 //      | AliFMDReconstructor |<>-----| AliFMDMultAlgorithm |
37 //      +---------------------+       +---------------------+
38 //                                               ^
39 //                                               |
40 //                                   +-----------+---------+
41 //                                   |                     |
42 //                         +-------------------+   +------------------+
43 //                         | AliFMDMultPoisson |   | AliFMDMultNaiive |
44 //                         +-------------------+   +------------------+
45 //
46 // AliFMDReconstructor acts as a manager class.  It contains a list of
47 // AliFMDMultAlgorithm objects.  The call graph looks something like 
48 //
49 //
50 //       +----------------------+            +----------------------+
51 //       | :AliFMDReconstructor |            | :AliFMDMultAlgorithm |
52 //       +----------------------+            +----------------------+
53 //                  |                                  |
54 //    Reconstruct  +-+                                 |
55 //    ------------>| |                         PreRun +-+
56 //                 | |------------------------------->| |   
57 //                 | |                                +-+
58 //                 | |-----+ (for each event)          |
59 //                 | |     | *ProcessEvent             |
60 //                 |+-+    |                           |
61 //                 || |<---+                 PreEvent +-+
62 //                 || |------------------------------>| |      
63 //                 || |                               +-+
64 //                 || |-----+                          |
65 //                 || |     | ProcessDigits            |
66 //                 ||+-+    |                          |
67 //                 ||| |<---+                          |
68 //                 ||| |         *ProcessDigit(digit) +-+
69 //                 ||| |----------------------------->| |
70 //                 ||| |                              +-+
71 //                 ||+-+                               |
72 //                 || |                     PostEvent +-+
73 //                 || |------------------------------>| |
74 //                 || |                               +-+
75 //                 |+-+                                |
76 //                 | |                        PostRun +-+
77 //                 | |------------------------------->| |
78 //                 | |                                +-+
79 //                 +-+                                 |
80 //                  |                                  |
81 //
82 //
83 // 
84 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
85 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
86 //
87 //
88 //____________________________________________________________________
89
90 #include <AliLog.h>                        // ALILOG_H
91 #include <AliRun.h>                        // ALIRUN_H
92 #include <AliRunLoader.h>                  // ALIRUNLOADER_H
93 #include <AliLoader.h>                     // ALILOADER_H
94 #include <AliHeader.h>                     // ALIHEADER_H
95 #include <AliRawReader.h>                  // ALIRAWREADER_H
96 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
97 #include "AliFMD.h"                        // ALIFMD_H
98 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
99 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
100 #include "AliFMDDetector.h"                // ALIFMDDETECTOR_H
101 #include "AliFMDRing.h"                    // ALIFMDRING_H
102 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
103 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
104 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
105 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
106 #include "AliFMDMultAlgorithm.h"           // ALIFMDMULTALGORITHM_H
107 #include "AliFMDMultPoisson.h"             // ALIFMDMULTPOISSON_H
108 #include "AliFMDMultNaiive.h"              // ALIFMDMULTNAIIVE_H
109
110 //____________________________________________________________________
111 ClassImp(AliFMDReconstructor)
112 #if 0
113   ; // This is here to keep Emacs for indenting the next line
114 #endif
115
116 //____________________________________________________________________
117 AliFMDReconstructor::AliFMDReconstructor() 
118   : AliReconstructor(),
119     fPedestal(0), 
120     fPedestalWidth(0),
121     fPedestalFactor(0)
122 {
123   // Make a new FMD reconstructor object - default CTOR.
124   AliFMDParameters* pars = AliFMDParameters::Instance();
125   fPedestal       = pars->GetPedestal();
126   fPedestalWidth  = pars->GetPedestalWidth();
127   fPedestalFactor = pars->GetPedestalFactor();
128   
129   fAlgorithms.Add(new AliFMDMultNaiive);
130   fAlgorithms.Add(new AliFMDMultPoisson);
131
132   TIter next(&fAlgorithms);
133   AliFMDMultAlgorithm* algorithm = 0;
134   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
135     algorithm->PreRun(0);
136 }
137   
138
139 //____________________________________________________________________
140 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
141   : AliReconstructor(),
142     fPedestal(0), 
143     fPedestalWidth(0),
144     fPedestalFactor(0)
145 {
146   // Copy constructor 
147   fPedestal       = other.fPedestal;
148   fPedestalWidth  = other.fPedestalWidth;
149   fPedestalFactor = other.fPedestalFactor;
150
151   
152   fAlgorithms.Delete();
153   TIter next(&(other.fAlgorithms));
154   AliFMDMultAlgorithm* algorithm = 0;
155   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
156     fAlgorithms.Add(algorithm);
157   fAlgorithms.SetOwner(kFALSE);
158 }
159   
160
161 //____________________________________________________________________
162 AliFMDReconstructor&
163 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
164 {
165   // Assignment operator
166   fPedestal       = other.fPedestal;
167   fPedestalWidth  = other.fPedestalWidth;
168   fPedestalFactor = other.fPedestalFactor;
169
170   fAlgorithms.Delete();
171   TIter next(&(other.fAlgorithms));
172   AliFMDMultAlgorithm* algorithm = 0;
173   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
174     fAlgorithms.Add(algorithm);
175   fAlgorithms.SetOwner(kFALSE);
176
177   return *this;
178 }
179
180 //____________________________________________________________________
181 AliFMDReconstructor::~AliFMDReconstructor() 
182 {
183   // Destructor 
184   fAlgorithms.Delete();
185 }
186
187 //____________________________________________________________________
188 void 
189 AliFMDReconstructor::Init(AliRunLoader* runLoader) 
190 {
191   // Initialize the reconstructor 
192   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
193   fCurrentVertex = 0;
194   if (!runLoader) { 
195     Warning("Init", "No run loader");
196     return;
197   }
198   AliHeader* header = runLoader->GetHeader();
199   if (!header) {
200     Warning("Init", "No header");
201     return;
202   }
203   AliGenEventHeader* eventHeader = header->GenEventHeader();
204   if (!eventHeader) {
205     Warning("Init", "no event header");
206     return;
207   }
208   TArrayF vtx;
209   eventHeader->PrimaryVertex(vtx);
210   fCurrentVertex = vtx[2];
211   AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f", 
212                    header->GetRun(), header->GetEvent(), fCurrentVertex));
213 }
214
215 //____________________________________________________________________
216 void 
217 AliFMDReconstructor::ConvertDigits(AliRawReader* reader, 
218                                    TTree* digitsTree) const
219 {
220   // Convert Raw digits to AliFMDDigit's in a tree 
221   AliFMDRawReader rawRead(reader, digitsTree);
222   // rawRead.SetSampleRate(fFMD->GetSampleRate());
223   rawRead.Exec();
224 }
225
226 //____________________________________________________________________
227 void 
228 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
229                                  TTree* clusterTree) const 
230 {
231   // Reconstruct event from digits in tree 
232   // Get the FMD branch holding the digits. 
233   // FIXME: The vertex may not be known yet, so we may have to move
234   // some of this to FillESD. 
235   AliDebug(1, "Reconstructing from digits in a tree");
236   
237   TBranch *digitBranch = digitsTree->GetBranch("FMD");
238   if (!digitBranch) {
239     Error("Exec", "No digit branch for the FMD found");
240     return;
241   }
242   TClonesArray* digits = new TClonesArray("AliFMDDigit");
243   digitBranch->SetAddress(&digits);
244
245   TIter next(&fAlgorithms);
246   AliFMDMultAlgorithm* algorithm = 0;
247   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
248     algorithm->PreEvent(clusterTree, fCurrentVertex);
249   digitBranch->GetEntry(0);
250   
251   ProcessDigits(digits);
252
253   next.Reset();
254   algorithm = 0;
255   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
256     algorithm->PostEvent();
257   clusterTree->Fill();
258   digits->Delete();
259   delete digits;
260 }
261  
262 //____________________________________________________________________
263 void
264 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
265 {
266   // For each digit, find the pseudo rapdity, azimuthal angle, and
267   // number of corrected ADC counts, and pass it on to the algorithms
268   // used. 
269   Int_t nDigits = digits->GetEntries();
270   AliDebug(1, Form("Got %d digits", nDigits));
271   for (Int_t i = 0; i < nDigits; i++) {
272     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
273     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
274     AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
275     if (!subDetector) { 
276       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
277       continue;
278     }
279     
280     AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
281     Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
282     if (!ring) {
283       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
284               digit->Ring());
285       break;
286     }
287     
288     Float_t  realZ    = fCurrentVertex + ringZ;
289     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
290                          / ring->GetNStrips() * (digit->Strip() + .5) 
291                          + ring->GetLowR());
292     Float_t  theta    = TMath::ATan2(stripR, realZ);
293     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
294                          * (digit->Sector() + .5));
295     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
296     UShort_t counts   = SubtractPedestal(digit);
297     
298     TIter next(&fAlgorithms);
299     AliFMDMultAlgorithm* algorithm = 0;
300     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
301       algorithm->ProcessDigit(digit, eta, phi, counts);
302   }
303 }
304       
305 //____________________________________________________________________
306 UShort_t
307 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
308 {
309   // Member function to subtract the pedestal from a digit
310   // This implementation does nothing, but a derived class could over
311   // load this to subtract a pedestal that was given in a database or
312   // something like that. 
313
314   Int_t counts = 0;
315   Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
316   if (digit->Count3() > 0)      counts = digit->Count3();
317   else if (digit->Count2() > 0) counts = digit->Count2();
318   else                          counts = digit->Count1();
319   counts = TMath::Max(Int_t(counts - ped), 0);
320   return  UShort_t(counts);
321 }
322
323 //____________________________________________________________________
324 void 
325 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
326                              TTree*  /* clusterTree */,
327                              AliESD* /* esd*/) const
328 {
329   // nothing to be done
330   // FIXME: The vertex may not be known when Reconstruct is executed,
331   // so we may have to move some of that member function here. 
332
333 }
334
335 //____________________________________________________________________
336 //
337 // EOF
338 //