Fixes, and extra debug
[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 #include "AliESD.h"                        // ALIESD_H
110
111 //____________________________________________________________________
112 ClassImp(AliFMDReconstructor)
113 #if 0
114   ; // This is here to keep Emacs for indenting the next line
115 #endif
116
117 //____________________________________________________________________
118 AliFMDReconstructor::AliFMDReconstructor() 
119   : AliReconstructor(),
120     fPedestal(0), 
121     fPedestalWidth(0),
122     fPedestalFactor(0)
123 {
124   // Make a new FMD reconstructor object - default CTOR.
125   AliFMDParameters* pars = AliFMDParameters::Instance();
126   fPedestal       = pars->GetPedestal();
127   fPedestalWidth  = pars->GetPedestalWidth();
128   fPedestalFactor = pars->GetPedestalFactor();
129   
130   fAlgorithms.Add(new AliFMDMultNaiive);
131   fAlgorithms.Add(new AliFMDMultPoisson);
132
133   TIter next(&fAlgorithms);
134   AliFMDMultAlgorithm* algorithm = 0;
135   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
136     algorithm->PreRun(0);
137 }
138   
139
140 //____________________________________________________________________
141 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
142   : AliReconstructor(),
143     fPedestal(0), 
144     fPedestalWidth(0),
145     fPedestalFactor(0)
146 {
147   // Copy constructor 
148   fPedestal       = other.fPedestal;
149   fPedestalWidth  = other.fPedestalWidth;
150   fPedestalFactor = other.fPedestalFactor;
151
152   
153   fAlgorithms.Delete();
154   TIter next(&(other.fAlgorithms));
155   AliFMDMultAlgorithm* algorithm = 0;
156   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
157     fAlgorithms.Add(algorithm);
158   fAlgorithms.SetOwner(kFALSE);
159 }
160   
161
162 //____________________________________________________________________
163 AliFMDReconstructor&
164 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
165 {
166   // Assignment operator
167   fPedestal       = other.fPedestal;
168   fPedestalWidth  = other.fPedestalWidth;
169   fPedestalFactor = other.fPedestalFactor;
170
171   fAlgorithms.Delete();
172   TIter next(&(other.fAlgorithms));
173   AliFMDMultAlgorithm* algorithm = 0;
174   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
175     fAlgorithms.Add(algorithm);
176   fAlgorithms.SetOwner(kFALSE);
177
178   return *this;
179 }
180
181 //____________________________________________________________________
182 AliFMDReconstructor::~AliFMDReconstructor() 
183 {
184   // Destructor 
185   fAlgorithms.Delete();
186 }
187
188 //____________________________________________________________________
189 void 
190 AliFMDReconstructor::Init(AliRunLoader* runLoader) 
191 {
192   // Initialize the reconstructor 
193   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
194   fCurrentVertex = 0;
195   if (!runLoader) { 
196     Warning("Init", "No run loader");
197     return;
198   }
199   AliHeader* header = runLoader->GetHeader();
200   if (!header) {
201     Warning("Init", "No header");
202     return;
203   }
204   AliGenEventHeader* eventHeader = header->GenEventHeader();
205   if (eventHeader) {
206     TArrayF vtx;
207     eventHeader->PrimaryVertex(vtx);
208     fCurrentVertex = vtx[2];
209     AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f", 
210                      header->GetRun(), header->GetEvent(), fCurrentVertex));
211     Warning("Init", "no generator event header");
212   }
213   else {
214     Warning("Init", "No generator event header - "
215             "perhaps we get the vertex from ESD?");
216   }
217   // Get the ESD tree 
218   SetESD(new AliESD);
219 }
220
221 //____________________________________________________________________
222 void 
223 AliFMDReconstructor::ConvertDigits(AliRawReader* reader, 
224                                    TTree* digitsTree) const
225 {
226   // Convert Raw digits to AliFMDDigit's in a tree 
227   AliFMDRawReader rawRead(reader, digitsTree);
228   // rawRead.SetSampleRate(fFMD->GetSampleRate());
229   rawRead.Exec();
230 }
231
232 //____________________________________________________________________
233 void 
234 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
235                                  TTree* clusterTree) const 
236 {
237   // Reconstruct event from digits in tree 
238   // Get the FMD branch holding the digits. 
239   // FIXME: The vertex may not be known yet, so we may have to move
240   // some of this to FillESD. 
241   AliDebug(1, "Reconstructing from digits in a tree");
242
243   if (fESD) {
244     const AliESDVertex* vertex = fESD->GetVertex();
245     // if (vertex) {
246     //   AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
247     //   fCurrentVertex = vertex->GetZv();
248     // }
249   }
250   
251   TBranch *digitBranch = digitsTree->GetBranch("FMD");
252   if (!digitBranch) {
253     Error("Exec", "No digit branch for the FMD found");
254     return;
255   }
256   TClonesArray* digits = new TClonesArray("AliFMDDigit");
257   digitBranch->SetAddress(&digits);
258
259   TIter next(&fAlgorithms);
260   AliFMDMultAlgorithm* algorithm = 0;
261   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
262     algorithm->PreEvent(clusterTree, fCurrentVertex);
263   digitBranch->GetEntry(0);
264   
265   ProcessDigits(digits);
266
267   next.Reset();
268   algorithm = 0;
269   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
270     algorithm->PostEvent();
271   clusterTree->Fill();
272   digits->Delete();
273   delete digits;
274 }
275  
276 //____________________________________________________________________
277 void
278 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
279 {
280   // For each digit, find the pseudo rapdity, azimuthal angle, and
281   // number of corrected ADC counts, and pass it on to the algorithms
282   // used. 
283   Int_t nDigits = digits->GetEntries();
284   AliDebug(1, Form("Got %d digits", nDigits));
285   for (Int_t i = 0; i < nDigits; i++) {
286     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
287     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
288     AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
289     if (!subDetector) { 
290       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
291       continue;
292     }
293     
294     AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
295     Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
296     if (!ring) {
297       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
298               digit->Ring());
299       break;
300     }
301     
302     Float_t  realZ    = fCurrentVertex + ringZ;
303     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
304                          / ring->GetNStrips() * (digit->Strip() + .5) 
305                          + ring->GetLowR());
306     Float_t  theta    = TMath::ATan2(stripR, realZ);
307     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
308                          * (digit->Sector() + .5));
309     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
310     UShort_t counts   = SubtractPedestal(digit);
311     
312     TIter next(&fAlgorithms);
313     AliFMDMultAlgorithm* algorithm = 0;
314     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
315       algorithm->ProcessDigit(digit, eta, phi, counts);
316   }
317 }
318       
319 //____________________________________________________________________
320 UShort_t
321 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
322 {
323   // Member function to subtract the pedestal from a digit
324   // This implementation does nothing, but a derived class could over
325   // load this to subtract a pedestal that was given in a database or
326   // something like that. 
327
328   Int_t counts = 0;
329   Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
330   if (digit->Count3() > 0)      counts = digit->Count3();
331   else if (digit->Count2() > 0) counts = digit->Count2();
332   else                          counts = digit->Count1();
333   counts = TMath::Max(Int_t(counts - ped), 0);
334   return  UShort_t(counts);
335 }
336
337 //____________________________________________________________________
338 void 
339 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
340                              TTree*  /* clusterTree */,
341                              AliESD* /* esd*/) const
342 {
343   // nothing to be done
344   // FIXME: The vertex may not be known when Reconstruct is executed,
345   // so we may have to move some of that member function here. 
346 #if 0
347   TClonesArray* multStrips  = 0;
348   TClonesArray* multRegions = 0;
349   TTree*        treeR  = fmdLoader->TreeR();
350   TBranch*      branchRegions = treeR->GetBranch("FMDPoisson");
351   TBranch*      branchStrips  = treeR->GetBranch("FMDNaiive");
352   branchRegions->SetAddress(&multRegions);
353   branchStrips->SetAddress(&multStrips);
354   
355   Int_t total = 0;
356   Int_t nEntries  = clusterTree->GetEntries();
357   for (Int_t entry = 0; entry < nEntries; entry++) {
358     AliDebug(5, Form("Entry # %d in cluster tree", entry));
359     treeR->GetEntry(entry);
360     
361     
362     Int_t nMults = multRegions->GetLast();
363     for (Int_t i = 0; i <= nMults; i++) {
364       AliFMDMultRegion* multR =
365         static_cast<AliFMDMultRegion*>(multRegions->UncheckedAt(i));
366       Int_t nParticles=multR->Particles();
367       if (i>=0 && i<=13)   hEtaPoissonI1->AddBinContent(i+1,nParticles);
368       if (i>=14 && i<=27 ) hEtaPoissonI2->AddBinContent(i-13,nParticles);
369       if (i>=28 && i<=33 );
370       if (i>=34 && i<=47 ) hEtaPoissonI3->AddBinContent(48-i,nParticles);
371       if (i>=48 && i<=53)  hEtaPoissonO3->AddBinContent(54-i,nParticles);
372     }
373   }
374 #endif   
375 }
376
377 //____________________________________________________________________
378 //
379 // EOF
380 //