Added docs and fixed a bug
[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 "AliFMDDigit.h"                   // ALIFMDDIGIT_H
99 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
100 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
101 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
102 #include "AliFMDMultAlgorithm.h"           // ALIFMDMULTALGORITHM_H
103 #include "AliFMDMultPoisson.h"             // ALIFMDMULTPOISSON_H
104 #include "AliFMDMultNaiive.h"              // ALIFMDMULTNAIIVE_H
105
106 //____________________________________________________________________
107 ClassImp(AliFMDReconstructor);
108
109 //____________________________________________________________________
110 AliFMDReconstructor::AliFMDReconstructor() 
111   : AliReconstructor(),
112     fPedestal(0), 
113     fPedestalWidth(0),
114     fPedestalFactor(0)
115 {
116   // Make a new FMD reconstructor object - default CTOR.
117   SetPedestal();
118
119   fFMDLoader = 0;
120   fRunLoader = 0;
121   fFMD       = 0;
122   fAlgorithms.Add(new AliFMDMultNaiive);
123   fAlgorithms.Add(new AliFMDMultPoisson);
124 }
125   
126
127 //____________________________________________________________________
128 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
129   : AliReconstructor(),
130     fPedestal(0), 
131     fPedestalWidth(0),
132     fPedestalFactor(0)
133 {
134   // Copy constructor 
135   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
136
137   fFMDLoader = other.fFMDLoader;
138   fRunLoader = other.fRunLoader;
139   fFMD       = other.fFMD;
140   
141   fAlgorithms.Delete();
142   TIter next(&(other.fAlgorithms));
143   AliFMDMultAlgorithm* algorithm = 0;
144   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
145     fAlgorithms.Add(algorithm);
146   fAlgorithms.SetOwner(kFALSE);
147 }
148   
149
150 //____________________________________________________________________
151 AliFMDReconstructor&
152 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
153 {
154   // Assignment operator
155   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
156
157   fFMDLoader = other.fFMDLoader;
158   fRunLoader = other.fRunLoader;
159   fFMD       = other.fFMD;
160
161   fAlgorithms.Delete();
162   TIter next(&(other.fAlgorithms));
163   AliFMDMultAlgorithm* algorithm = 0;
164   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
165     fAlgorithms.Add(algorithm);
166   fAlgorithms.SetOwner(kFALSE);
167
168   return *this;
169 }
170
171 //____________________________________________________________________
172 AliFMDReconstructor::~AliFMDReconstructor() 
173 {
174   // Destructor 
175   fAlgorithms.Delete();
176 }
177   
178 //____________________________________________________________________
179 void 
180 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor) 
181 {
182   // Set the pedestal, and pedestal width 
183   fPedestal       = mean;
184   fPedestalWidth  = width;
185   fPedestalFactor = factor;
186 }
187
188 //____________________________________________________________________
189 void 
190 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader, 
191                                  AliRawReader* rawReader) const
192
193   // Collects all digits in the same active volume into number of
194   // particles
195   //
196   // Reconstruct number of particles in given group of pads for given
197   // FMDvolume determined by numberOfVolume,
198   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
199   // numberOgMaxRing 
200   //
201   // The reconstruction method is choosen based on the number of empty
202   // strips. 
203   if (!runLoader) {
204     Error("Exec","Run Loader loader is NULL - Session not opened");
205     return;
206   }
207   fRunLoader = runLoader;
208   fFMDLoader = runLoader->GetLoader("FMDLoader");
209   if (!fFMDLoader) 
210     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
211           "in specified event");
212
213   // Get the AliRun object
214   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
215   
216   // Get the AliFMD object
217   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
218   if (!fFMD) {
219     AliError("Can not get FMD from gAlice");
220     return;
221   }
222   fFMDLoader->LoadRecPoints("RECREATE");
223
224   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
225
226   TIter next(&fAlgorithms);
227   AliFMDMultAlgorithm* algorithm = 0;
228   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
229     algorithm->PreRun(fFMD);
230
231   if (rawReader) {
232     Int_t event = 0;
233     while (rawReader->NextEvent()) {
234       ProcessEvent(event, rawReader);
235       event++;
236     }
237   }
238   else {
239     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries()); 
240     for(Int_t event = 0; event < nEvents; event++) 
241       ProcessEvent(event, 0);
242   }
243
244   next.Reset();
245   algorithm = 0;
246   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
247     algorithm->PostRun();
248
249   fFMDLoader->UnloadRecPoints();
250   fFMDLoader = 0;
251   fRunLoader = 0;
252   fFMD       = 0;
253 }
254
255 //____________________________________________________________________
256 void 
257 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
258
259   // Collects all digits in the same active volume into number of
260   // particles
261   //
262   // Reconstruct number of particles in given group of pads for given
263   // FMDvolume determined by numberOfVolume,
264   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
265   // numberOgMaxRing 
266   //
267   // The reconstruction method is choosen based on the number of empty
268   // strips. 
269   Reconstruct(runLoader, 0);
270 }
271
272
273 //____________________________________________________________________
274 void 
275 AliFMDReconstructor::ProcessEvent(Int_t event, 
276                                   AliRawReader* reader) const
277 {
278   // Process one event read from either a clones array or from a a raw
279   // data reader. 
280   fRunLoader->GetEvent(event) ;
281   //event z-vertex for correction eta-rad dependence      
282   AliHeader *header            = fRunLoader->GetHeader();
283   if (!header) Warning("ProcessEvent", "no AliHeader found!");
284   AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
285
286   // Get the Z--coordinate from the event header 
287   TArrayF o(3); 
288   if (genHeader) genHeader->PrimaryVertex(o);
289   else           Warning("ProcessEvent", "No GenEventHeader Found");
290   fCurrentVertex = o.At(2);
291
292   // If the recontruction tree isn't loaded, load it
293   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
294   
295   // Load or recreate the digits 
296   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
297     if (!reader) {
298       Error("Exec","Error occured while loading digits. Exiting.");
299       return;
300     }
301   }
302   // Get the digits tree 
303   TTree* digitTree = fFMDLoader->TreeD();
304   if (!digitTree) { 
305     if (!reader) {
306       Error("Exec","Can not get Tree with Digits. "
307             "Nothing to reconstruct - Exiting");
308       return;
309     }
310     fFMDLoader->MakeTree("D");
311     digitTree = fFMDLoader->TreeD();
312     
313   }
314   // Get the FMD branch holding the digits. 
315   TBranch *digitBranch = digitTree->GetBranch("FMD");
316   TClonesArray* digits = fFMD->Digits();
317   if (!digitBranch) {
318     if (!reader) {
319       Error("Exec", "No digit branch for the FMD found");
320       return;
321     }
322     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
323   }
324   if (!reader) digitBranch->SetAddress(&digits);
325
326   if  (reader) {
327     AliFMDRawReader rawRead(fFMD, reader);
328     AliDebug(10, Form("Making raw reader with sample rate: %d",
329                       fFMD->GetSampleRate()));
330     rawRead.SetSampleRate(fFMD->GetSampleRate());
331     rawRead.Exec();
332   }
333   else {
334     // Read the ADC values from a clones array. 
335     AliDebug(10, "Reading ADCs from Digits array");
336     // read Digits, and reconstruct the particles
337     if (!fFMDLoader->TreeD()->GetEvent(0)) return;
338   }
339   
340   TIter next(&fAlgorithms);
341   AliFMDMultAlgorithm* algorithm = 0;
342   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
343     algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);
344
345   ProcessDigits(digits);
346
347   next.Reset();
348   algorithm = 0;
349   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
350     algorithm->PostEvent();
351   
352   if (reader) {
353     digitTree->Fill();
354     fFMDLoader->WriteDigits("OVERWRITE");
355   }
356   fFMDLoader->UnloadDigits();
357   fFMDLoader->TreeR()->Reset();
358   fFMDLoader->TreeR()->Fill(); 
359   fFMDLoader->WriteRecPoints("OVERWRITE");
360 }
361
362 //____________________________________________________________________
363 UShort_t
364 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
365 {
366   // Member function to subtract the pedestal from a digit
367   // This implementation does nothing, but a derived class could over
368   // load this to subtract a pedestal that was given in a database or
369   // something like that. 
370
371   Int_t counts = 0;
372   Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
373   if (digit->Count3() > 0)      counts = digit->Count3();
374   else if (digit->Count2() > 0) counts = digit->Count2();
375   else                          counts = digit->Count1();
376   counts = TMath::Max(Int_t(counts - ped), 0);
377   return  UShort_t(counts);
378 }
379
380 //____________________________________________________________________
381 void
382 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
383 {
384   Int_t nDigits = digits->GetEntries();
385   for (Int_t i = 0; i < nDigits; i++) {
386     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
387     AliFMDSubDetector* subDetector = 0;
388     switch (digit->Detector()) {
389     case 1: subDetector = fFMD->GetFMD1(); break;
390     case 2: subDetector = fFMD->GetFMD2(); break;
391     case 3: subDetector = fFMD->GetFMD3(); break;
392     }
393     if (!subDetector) { 
394       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
395       continue;
396     }
397     
398     AliFMDRing* ring  = 0;
399     Float_t     ringZ = 0;
400     switch(digit->Ring()) {
401     case 'i':
402     case 'I':  
403       ring  = subDetector->GetInner(); 
404       ringZ = subDetector->GetInnerZ();
405       break;
406     case 'o':
407     case 'O':  
408       ring  = subDetector->GetOuter(); 
409       ringZ = subDetector->GetOuterZ();
410       break;
411     }
412     if (!ring) {
413       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
414               digit->Ring());
415       break;
416     }
417     
418     Float_t  realZ    = fCurrentVertex + ringZ;
419     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
420                          / ring->GetNStrips() * (digit->Strip() + .5) 
421                          + ring->GetLowR());
422     Float_t  theta    = TMath::ATan2(stripR, realZ);
423     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
424                          * (digit->Sector() + .5));
425     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
426     UShort_t counts   = SubtractPedestal(digit);
427     
428     TIter next(&fAlgorithms);
429     AliFMDMultAlgorithm* algorithm = 0;
430     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
431       algorithm->ProcessDigit(digit, eta, phi, counts);
432   }
433 }
434       
435  
436 //____________________________________________________________________
437 void 
438 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
439                              AliESD* /*esd*/) const
440 {
441   // nothing to be done
442
443 }
444
445 //____________________________________________________________________
446 //
447 // EOF
448 //