49ea78a70a024b58a020eb57917e67ca3a97babb
[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 ReconstParticles (reconstructed
21 // particles) out of Digits
22 //
23 //
24 // This class reads either digits from a TClonesArray or raw data from
25 // a DDL file (or similar), and stores the read ADC counts in an
26 // internal cache (fAdcs). 
27 //
28 // From the cached values it then calculates the number of particles
29 // that hit a region of the FMDs, as specified by the user. 
30 //
31 // The reconstruction can be done in two ways: Either via counting the
32 // number of empty strips (Poisson method), or by converting the ADC
33 // signal to an energy deposition, and then dividing by the typical
34 // energy loss of a particle.
35 // 
36 // Currently, this class only reads the digits from a TClonesArray,
37 // and the Poission method for reconstruction. 
38 //
39 // 
40 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
41 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
42 //
43 //
44 //____________________________________________________________________
45
46 #include <AliLog.h>                        // ALILOG_H
47 #include <AliRun.h>                        // ALIRUN_H
48 #include <AliRunLoader.h>                  // ALIRUNLOADER_H
49 #include <AliLoader.h>                     // ALILOADER_H
50 #include <AliHeader.h>                     // ALIHEADER_H
51 #include <AliRawReader.h>                  // ALIRAWREADER_H
52 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
53 #include "AliFMD.h"                        // ALIFMD_H
54 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
55 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
56 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
57 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
58 #include "AliFMDMultAlgorithm.h"           // ALIFMDMULTALGORITHM_H
59 #include "AliFMDMultPoisson.h"             // ALIFMDMULTPOISSON_H
60 #include "AliFMDMultNaiive.h"              // ALIFMDMULTNAIIVE_H
61
62 //____________________________________________________________________
63 ClassImp(AliFMDReconstructor);
64
65 //____________________________________________________________________
66 AliFMDReconstructor::AliFMDReconstructor() 
67   : AliReconstructor(),
68     fPedestal(0), 
69     fPedestalWidth(0),
70     fPedestalFactor(0)
71 {
72   // Make a new FMD reconstructor object - default CTOR.
73   SetPedestal();
74
75   fFMDLoader = 0;
76   fRunLoader = 0;
77   fFMD       = 0;
78   fAlgorithms.Add(new AliFMDMultNaiive);
79   fAlgorithms.Add(new AliFMDMultPoisson);
80 }
81   
82
83 //____________________________________________________________________
84 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
85   : AliReconstructor(),
86     fPedestal(0), 
87     fPedestalWidth(0),
88     fPedestalFactor(0)
89 {
90   // Copy constructor 
91   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
92
93   fFMDLoader = other.fFMDLoader;
94   fRunLoader = other.fRunLoader;
95   fFMD       = other.fFMD;
96   
97   fAlgorithms.Delete();
98   TIter next(&(other.fAlgorithms));
99   AliFMDMultAlgorithm* algorithm = 0;
100   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
101     fAlgorithms.Add(algorithm);
102   fAlgorithms.SetOwner(kFALSE);
103 }
104   
105
106 //____________________________________________________________________
107 AliFMDReconstructor&
108 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
109 {
110   // Assignment operator
111   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
112
113   fFMDLoader = other.fFMDLoader;
114   fRunLoader = other.fRunLoader;
115   fFMD       = other.fFMD;
116
117   fAlgorithms.Delete();
118   TIter next(&(other.fAlgorithms));
119   AliFMDMultAlgorithm* algorithm = 0;
120   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
121     fAlgorithms.Add(algorithm);
122   fAlgorithms.SetOwner(kFALSE);
123
124   return *this;
125 }
126
127 //____________________________________________________________________
128 AliFMDReconstructor::~AliFMDReconstructor() 
129 {
130   // Destructor 
131   fAlgorithms.Delete();
132 }
133   
134 //____________________________________________________________________
135 void 
136 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor) 
137 {
138   // Set the pedestal, and pedestal width 
139   fPedestal       = mean;
140   fPedestalWidth  = width;
141   fPedestalFactor = factor;
142 }
143
144 //____________________________________________________________________
145 void 
146 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader, 
147                                  AliRawReader* rawReader) const
148
149   // Collects all digits in the same active volume into number of
150   // particles
151   //
152   // Reconstruct number of particles in given group of pads for given
153   // FMDvolume determined by numberOfVolume,
154   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
155   // numberOgMaxRing 
156   //
157   // The reconstruction method is choosen based on the number of empty
158   // strips. 
159   if (!runLoader) {
160     Error("Exec","Run Loader loader is NULL - Session not opened");
161     return;
162   }
163   fRunLoader = runLoader;
164   fFMDLoader = runLoader->GetLoader("FMDLoader");
165   if (!fFMDLoader) 
166     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
167           "in specified event");
168
169   // Get the AliRun object
170   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
171   
172   // Get the AliFMD object
173   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
174   if (!fFMD) {
175     AliError("Can not get FMD from gAlice");
176     return;
177   }
178   fFMDLoader->LoadRecPoints("RECREATE");
179
180   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
181
182   TIter next(&fAlgorithms);
183   AliFMDMultAlgorithm* algorithm = 0;
184   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
185     algorithm->PreRun(fFMD);
186
187   if (rawReader) {
188     Int_t event = 0;
189     while (rawReader->NextEvent()) {
190       ProcessEvent(event, rawReader);
191       event++;
192     }
193   }
194   else {
195     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries()); 
196     for(Int_t event = 0; event < nEvents; event++) 
197       ProcessEvent(event, 0);
198   }
199
200   next.Reset();
201   algorithm = 0;
202   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
203     algorithm->PostRun();
204
205   fFMDLoader->UnloadRecPoints();
206   fFMDLoader = 0;
207   fRunLoader = 0;
208   fFMD       = 0;
209 }
210
211 //____________________________________________________________________
212 void 
213 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
214
215   // Collects all digits in the same active volume into number of
216   // particles
217   //
218   // Reconstruct number of particles in given group of pads for given
219   // FMDvolume determined by numberOfVolume,
220   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
221   // numberOgMaxRing 
222   //
223   // The reconstruction method is choosen based on the number of empty
224   // strips. 
225   Reconstruct(runLoader, 0);
226 }
227
228
229 //____________________________________________________________________
230 void 
231 AliFMDReconstructor::ProcessEvent(Int_t event, 
232                                   AliRawReader* reader) const
233 {
234   // Process one event read from either a clones array or from a a raw
235   // data reader. 
236   fRunLoader->GetEvent(event) ;
237   //event z-vertex for correction eta-rad dependence      
238   AliHeader *header            = fRunLoader->GetHeader();
239   if (!header) Warning("ProcessEvent", "no AliHeader found!");
240   AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
241
242   // Get the Z--coordinate from the event header 
243   TArrayF o(3); 
244   if (genHeader) genHeader->PrimaryVertex(o);
245   else           Warning("ProcessEvent", "No GenEventHeader Found");
246   fCurrentVertex = o.At(2);
247
248   // If the recontruction tree isn't loaded, load it
249   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
250   
251   // Load or recreate the digits 
252   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
253     if (!reader) {
254       Error("Exec","Error occured while loading digits. Exiting.");
255       return;
256     }
257   }
258   // Get the digits tree 
259   TTree* digitTree = fFMDLoader->TreeD();
260   if (!digitTree) { 
261     if (!reader) {
262       Error("Exec","Can not get Tree with Digits. "
263             "Nothing to reconstruct - Exiting");
264       return;
265     }
266     fFMDLoader->MakeTree("D");
267     digitTree = fFMDLoader->TreeD();
268     
269   }
270   // Get the FMD branch holding the digits. 
271   TBranch *digitBranch = digitTree->GetBranch("FMD");
272   TClonesArray* digits = fFMD->Digits();
273   if (!digitBranch) {
274     if (!reader) {
275       Error("Exec", "No digit branch for the FMD found");
276       return;
277     }
278     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
279   }
280   if (!reader) digitBranch->SetAddress(&digits);
281
282   if  (reader) {
283     AliFMDRawReader rawRead(fFMD, reader);
284     // rawRead->SetSampleRate(fSampleRate);
285     rawRead.Exec();
286   }
287   else {
288     // Read the ADC values from a clones array. 
289     AliDebug(10, "Reading ADCs from Digits array");
290     // read Digits, and reconstruct the particles
291     if (!fFMDLoader->TreeD()->GetEvent(0)) return;
292   }
293   
294   TIter next(&fAlgorithms);
295   AliFMDMultAlgorithm* algorithm = 0;
296   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
297     algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);
298
299   ProcessDigits(digits);
300
301   next.Reset();
302   algorithm = 0;
303   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
304     algorithm->PostEvent();
305   
306   if (reader) {
307     digitTree->Fill();
308     fFMDLoader->WriteDigits("OVERWRITE");
309   }
310   fFMDLoader->UnloadDigits();
311   fFMDLoader->TreeR()->Reset();
312   fFMDLoader->TreeR()->Fill(); 
313   fFMDLoader->WriteRecPoints("OVERWRITE");
314 }
315
316 //____________________________________________________________________
317 UShort_t
318 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
319 {
320   // Member function to subtract the pedestal from a digit
321   // This implementation does nothing, but a derived class could over
322   // load this to subtract a pedestal that was given in a database or
323   // something like that. 
324
325   Int_t counts = 0;
326   Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
327   if (digit->Count3() >= 0)      counts = digit->Count3();
328   else if (digit->Count2() >= 0) counts = digit->Count2();
329   else                           counts = digit->Count2();
330   counts = TMath::Max(Int_t(counts - ped), 0);
331   return  UShort_t(counts);
332 }
333
334 //____________________________________________________________________
335 void
336 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
337 {
338   Int_t nDigits = digits->GetEntries();
339   for (Int_t i = 0; i < nDigits; i++) {
340     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
341     AliFMDSubDetector* subDetector = 0;
342     switch (digit->Detector()) {
343     case 1: subDetector = fFMD->GetFMD1(); break;
344     case 2: subDetector = fFMD->GetFMD2(); break;
345     case 3: subDetector = fFMD->GetFMD3(); break;
346     }
347     if (!subDetector) { 
348       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
349       continue;
350     }
351     
352     AliFMDRing* ring  = 0;
353     Float_t     ringZ = 0;
354     switch(digit->Ring()) {
355     case 'i':
356     case 'I':  
357       ring  = subDetector->GetInner(); 
358       ringZ = subDetector->GetInnerZ();
359       break;
360     case 'o':
361     case 'O':  
362       ring  = subDetector->GetOuter(); 
363       ringZ = subDetector->GetOuterZ();
364       break;
365     }
366     if (!ring) {
367       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
368               digit->Ring());
369       break;
370     }
371     
372     Float_t  realZ    = fCurrentVertex + ringZ;
373     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
374                          / ring->GetNStrips() * (digit->Strip() + .5) 
375                          + ring->GetLowR());
376     Float_t  theta    = TMath::ATan2(stripR, realZ);
377     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
378                          * (digit->Sector() + .5));
379     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
380     UShort_t counts   = SubtractPedestal(digit);
381     
382     TIter next(&fAlgorithms);
383     AliFMDMultAlgorithm* algorithm = 0;
384     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
385       algorithm->ProcessDigit(digit, eta, phi, counts);
386   }
387 }
388       
389  
390 //____________________________________________________________________
391 void 
392 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/, 
393                              AliESD* /*esd*/) const
394 {
395   // nothing to be done
396
397 }
398
399 //____________________________________________________________________
400 //
401 // EOF
402 //