The patch fixes
[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 /* $Id$ */
16 /** @file    AliFMDReconstructor.cxx
17     @author  Christian Holm Christensen <cholm@nbi.dk>
18     @date    Mon Mar 27 12:47:09 2006
19     @brief   FMD reconstruction 
20 */
21 //____________________________________________________________________
22 //
23 // This is a class that constructs AliFMDRecPoint objects from of Digits
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).   The rec-points are made via the naiive
27 // method. 
28 //
29 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
30 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
31 //
32 //
33 //____________________________________________________________________
34
35 #include <AliLog.h>                        // ALILOG_H
36 // #include <AliRun.h>                        // ALIRUN_H
37 #include <AliRunLoader.h>                  // ALIRUNLOADER_H
38 #include <AliHeader.h>                     // ALIHEADER_H
39 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
40 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
41 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
42 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
43 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
44 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
45 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
46 #include "AliESD.h"                        // ALIESD_H
47 #include <AliESDFMD.h>                     // ALIESDFMD_H
48 class AliRawReader;
49
50 //____________________________________________________________________
51 ClassImp(AliFMDReconstructor)
52 #if 0
53   ; // This is here to keep Emacs for indenting the next line
54 #endif
55
56 //____________________________________________________________________
57 AliFMDReconstructor::AliFMDReconstructor() 
58   : AliReconstructor(),
59     fMult(0x0), 
60     fNMult(0),
61     fTreeR(0x0),
62     fCurrentVertex(0),
63     fESDObj(0x0),
64     fESD(0x0)
65 {
66   // Make a new FMD reconstructor object - default CTOR.  
67 }
68   
69
70 //____________________________________________________________________
71 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
72   : AliReconstructor(), 
73     fMult(other.fMult),
74     fNMult(other.fNMult),
75     fTreeR(other.fTreeR),
76     fCurrentVertex(other.fCurrentVertex),
77     fESDObj(other.fESDObj),
78     fESD(other.fESD)
79 {
80   // Copy constructor 
81 }
82   
83
84 //____________________________________________________________________
85 AliFMDReconstructor&
86 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
87 {
88   // Assignment operator
89   fMult   = other.fMult;
90   fNMult = other.fNMult;
91   fTreeR = other.fTreeR;
92   fCurrentVertex = other.fCurrentVertex;
93   fESDObj = other.fESDObj;
94   fESD = other.fESD;
95   return *this;
96 }
97
98 //____________________________________________________________________
99 AliFMDReconstructor::~AliFMDReconstructor() 
100 {
101   // Destructor 
102   if (fMult)   fMult->Delete();
103   if (fMult)   delete fMult;
104   if (fESDObj) delete fESDObj;
105 }
106
107 //____________________________________________________________________
108 void 
109 AliFMDReconstructor::Init(AliRunLoader* runLoader) 
110 {
111   // Initialize the reconstructor 
112   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
113   // Initialize the geometry 
114   AliFMDGeometry* geom = AliFMDGeometry::Instance();
115   geom->Init();
116   geom->InitTransformations();
117
118   // Initialize the parameters
119   AliFMDParameters* param = AliFMDParameters::Instance();
120   param->Init();
121   
122   // Current vertex position
123   fCurrentVertex = 0;
124   // Create array of reconstructed strip multiplicities 
125   fMult = new TClonesArray("AliFMDRecPoint", 51200);
126   // Create ESD output object 
127   fESDObj = new AliESDFMD;
128   
129   // Check that we have a run loader
130   if (!runLoader) { 
131     Warning("Init", "No run loader");
132     return;
133   }
134
135   // Check if we can get the header tree 
136   AliHeader* header = runLoader->GetHeader();
137   if (!header) {
138     Warning("Init", "No header");
139     return;
140   }
141
142   // Check if we can get a simulation header 
143   AliGenEventHeader* eventHeader = header->GenEventHeader();
144   if (eventHeader) {
145     TArrayF vtx;
146     eventHeader->PrimaryVertex(vtx);
147     fCurrentVertex = vtx[2];
148     AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f", 
149                      header->GetRun(), header->GetEvent(), fCurrentVertex));
150     Warning("Init", "no generator event header");
151   }
152   else {
153     Warning("Init", "No generator event header - "
154             "perhaps we get the vertex from ESD?");
155   }
156 }
157
158 //____________________________________________________________________
159 void 
160 AliFMDReconstructor::ConvertDigits(AliRawReader* reader, 
161                                    TTree* digitsTree) const
162 {
163   // Convert Raw digits to AliFMDDigit's in a tree 
164   AliDebug(1, "Reading raw data into digits tree");
165   AliFMDRawReader rawRead(reader, digitsTree);
166   // rawRead.SetSampleRate(fFMD->GetSampleRate());
167   rawRead.Exec();
168 }
169
170 //____________________________________________________________________
171 void 
172 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
173                                  TTree* clusterTree) const 
174 {
175   // Reconstruct event from digits in tree 
176   // Get the FMD branch holding the digits. 
177   // FIXME: The vertex may not be known yet, so we may have to move
178   // some of this to FillESD. 
179   AliDebug(1, "Reconstructing from digits in a tree");
180 #if 1
181   if (fESD) {
182     const AliESDVertex* vertex = fESD->GetVertex();
183     if (vertex) {
184       AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
185       fCurrentVertex = vertex->GetZv();
186     }
187   }
188 #endif  
189   TBranch *digitBranch = digitsTree->GetBranch("FMD");
190   if (!digitBranch) {
191     Error("Exec", "No digit branch for the FMD found");
192     return;
193   }
194   TClonesArray* digits = new TClonesArray("AliFMDDigit");
195   digitBranch->SetAddress(&digits);
196
197   if (fMult)   fMult->Clear();
198   if (fESDObj) fESDObj->Clear();
199   
200   fNMult = 0;
201   fTreeR = clusterTree;
202   fTreeR->Branch("FMD", &fMult);
203   
204   AliDebug(5, "Getting entry 0 from digit branch");
205   digitBranch->GetEntry(0);
206   
207   AliDebug(5, "Processing digits");
208   ProcessDigits(digits);
209
210   Int_t written = clusterTree->Fill();
211   AliDebug(10, Form("Filled %d bytes into cluster tree", written));
212   digits->Delete();
213   delete digits;
214 }
215  
216
217 //____________________________________________________________________
218 void
219 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
220 {
221   // For each digit, find the pseudo rapdity, azimuthal angle, and
222   // number of corrected ADC counts, and pass it on to the algorithms
223   // used. 
224   Int_t nDigits = digits->GetEntries();
225   AliDebug(1, Form("Got %d digits", nDigits));
226   for (Int_t i = 0; i < nDigits; i++) {
227     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
228     AliFMDParameters* param  = AliFMDParameters::Instance();
229     // Check that the strip is not marked as dead 
230     if (param->IsDead(digit->Detector(), digit->Ring(), 
231                       digit->Sector(), digit->Strip())) {
232       AliDebug(10, Form("FMD%d%c[%2d,%3d] is dead", digit->Detector(), 
233                         digit->Ring(), digit->Sector(), digit->Strip()));
234       continue;
235     }
236
237     // digit->Print();
238     // Get eta and phi 
239     Float_t eta, phi;
240     PhysicalCoordinates(digit, eta, phi);
241     
242     // Substract pedestal. 
243     UShort_t counts   = SubtractPedestal(digit);
244     
245     // Gain match digits. 
246     Double_t edep     = Adc2Energy(digit, eta, counts);
247     
248     // Make rough multiplicity 
249     Double_t mult     = Energy2Multiplicity(digit, edep);
250     
251     AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
252                       "ADC: %d, Counts: %d, Energy: %f, Mult: %f", 
253                       digit->Detector(), digit->Ring(), digit->Sector(),
254                       digit->Strip(), digit->Counts(), counts, edep, mult));
255     
256     // Create a `RecPoint' on the output branch. 
257     AliFMDRecPoint* m = 
258       new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(), 
259                                             digit->Ring(), 
260                                             digit->Sector(),
261                                             digit->Strip(),
262                                             eta, phi, 
263                                             edep, mult);
264     (void)m; // Suppress warnings about unused variables. 
265     fNMult++;
266
267     fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(), 
268                              digit->Sector(),  digit->Strip(), mult);
269     fESDObj->SetEta(digit->Detector(), digit->Ring(), 
270                     digit->Sector(),  digit->Strip(), eta);
271   }
272 }
273
274 //____________________________________________________________________
275 UShort_t
276 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
277 {
278   // Member function to subtract the pedestal from a digit
279   // This implementation does nothing, but a derived class could over
280   // load this to subtract a pedestal that was given in a database or
281   // something like that. 
282
283   Int_t             counts = 0;
284   AliFMDParameters* param  = AliFMDParameters::Instance();
285   Float_t           pedM   = param->GetPedestal(digit->Detector(), 
286                                                 digit->Ring(), 
287                                                 digit->Sector(), 
288                                                 digit->Strip());
289   AliDebug(15, Form("Subtracting pedestal %f from signal %d", 
290                    pedM, digit->Counts()));
291   if (digit->Count3() > 0)      counts = digit->Count3();
292   else if (digit->Count2() > 0) counts = digit->Count2();
293   else                          counts = digit->Count1();
294   counts = TMath::Max(Int_t(counts - pedM), 0);
295   if (counts > 0) AliDebug(15, "Got a hit strip");
296   
297   return  UShort_t(counts);
298 }
299
300 //____________________________________________________________________
301 Float_t
302 AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit, 
303                                 Float_t      /* eta */, 
304                                 UShort_t     count) const
305 {
306   // Converts number of ADC counts to energy deposited. 
307   // Note, that this member function can be overloaded by derived
308   // classes to do strip-specific look-ups in databases or the like,
309   // to find the proper gain for a strip. 
310   // 
311   // In this simple version, we calculate the energy deposited as 
312   // 
313   //    EnergyDeposited = cos(theta) * gain * count
314   // 
315   // where 
316   // 
317   //           Pre_amp_MIP_Range
318   //    gain = ----------------- * Energy_deposited_per_MIP
319   //           ADC_channel_size    
320   // 
321   // is constant and the same for all strips.
322
323   // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
324   // Double_t edep  = TMath::Abs(TMath::Cos(theta)) * fGain * count;
325   AliFMDParameters* param = AliFMDParameters::Instance();
326   Float_t           gain  = param->GetPulseGain(digit->Detector(), 
327                                                 digit->Ring(), 
328                                                 digit->Sector(), 
329                                                 digit->Strip());
330   Double_t          edep  = count * gain;
331   AliDebug(15, Form("Converting counts %d to energy via factor %f", 
332                     count, gain));
333   return edep;
334 }
335
336 //____________________________________________________________________
337 Float_t
338 AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */, 
339                                          Float_t      edep) const
340 {
341   // Converts an energy signal to number of particles. 
342   // Note, that this member function can be overloaded by derived
343   // classes to do strip-specific look-ups in databases or the like,
344   // to find the proper gain for a strip. 
345   // 
346   // In this simple version, we calculate the multiplicity as 
347   // 
348   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
349   // 
350   // where 
351   //
352   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
353   // 
354   // is constant and the same for all strips 
355   AliFMDParameters* param   = AliFMDParameters::Instance();
356   Double_t          edepMIP = param->GetEdepMip();
357   Float_t           mult    = edep / edepMIP;
358   if (edep > 0) 
359     AliDebug(15, Form("Translating energy %f to multiplicity via "
360                      "divider %f->%f", edep, edepMIP, mult));
361   return mult;
362 }
363
364 //____________________________________________________________________
365 void
366 AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit, 
367                                          Float_t& eta, 
368                                          Float_t& phi) const
369 {
370   // Get the eta and phi of a digit 
371   // 
372   // Get geometry. 
373   AliFMDGeometry* geom = AliFMDGeometry::Instance();
374   Double_t x, y, z, r, theta;
375   geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(), 
376                     digit->Strip(), x, y, z);
377   // Correct for vertex offset. 
378   z     += fCurrentVertex;
379   phi   =  TMath::ATan2(y, x);
380   r     =  TMath::Sqrt(y * y + x * x);
381   theta =  TMath::ATan2(r, z);
382   eta   = -TMath::Log(TMath::Tan(theta / 2));
383 }
384
385       
386
387 //____________________________________________________________________
388 void 
389 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
390                              TTree*  /* clusterTree */,
391                              AliESD* esd) const
392 {
393   // nothing to be done
394   // FIXME: The vertex may not be known when Reconstruct is executed,
395   // so we may have to move some of that member function here. 
396   AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
397   // fESDObj->Print();
398
399   if (esd) { 
400     AliDebug(1, Form("Writing FMD data to ESD tree"));
401     esd->SetFMDData(fESDObj);
402   }
403 }
404
405
406 //____________________________________________________________________
407 void 
408 AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const 
409 {
410   // Cannot be used.  See member function with same name but with 2
411   // TTree arguments.   Make sure you do local reconstrucion 
412   AliDebug(1, Form("Calling FillESD with loader and tree"));
413   AliError("MayNotUse");
414 }
415 //____________________________________________________________________
416 void 
417 AliFMDReconstructor::Reconstruct(AliRunLoader*) const 
418 {
419   // Cannot be used.  See member function with same name but with 2
420   // TTree arguments.   Make sure you do local reconstrucion 
421   AliDebug(1, Form("Calling FillESD with loader"));
422   AliError("MayNotUse");
423 }
424 //____________________________________________________________________
425 void 
426 AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const 
427 {
428   // Cannot be used.  See member function with same name but with 2
429   // TTree arguments.   Make sure you do local reconstrucion 
430   AliDebug(1, Form("Calling FillESD with loader and raw reader"));
431   AliError("MayNotUse");
432 }
433 //____________________________________________________________________
434 void 
435 AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const 
436 {
437   // Cannot be used.  See member function with same name but with 2
438   // TTree arguments.   Make sure you do local reconstrucion 
439   AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
440   AliError("MayNotUse");
441 }
442 //____________________________________________________________________
443 void 
444 AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
445 {
446   // Cannot be used.  See member function with same name but with 2
447   // TTree arguments.   Make sure you do local reconstrucion 
448   AliDebug(1, Form("Calling FillESD with loader and ESD"));
449   AliError("MayNotUse");
450 }
451 //____________________________________________________________________
452 void 
453 AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const 
454 {
455   // Cannot be used.  See member function with same name but with 2
456   // TTree arguments.   Make sure you do local reconstrucion 
457   AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
458   AliError("MayNotUse");
459 }
460
461 //____________________________________________________________________
462 //
463 // EOF
464 //