ea8e2cd4d3da80ca6663aaad5b9d7485cb9cebaa
[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())) continue;
232
233     // digit->Print();
234     // Get eta and phi 
235     Float_t eta, phi;
236     PhysicalCoordinates(digit, eta, phi);
237     
238     // Substract pedestal. 
239     UShort_t counts   = SubtractPedestal(digit);
240     
241     // Gain match digits. 
242     Double_t edep     = Adc2Energy(digit, eta, counts);
243     
244     // Make rough multiplicity 
245     Double_t mult     = Energy2Multiplicity(digit, edep);
246     
247     AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
248                       "ADC: %d, Counts: %d, Energy: %f, Mult: %f", 
249                       digit->Detector(), digit->Ring(), digit->Sector(),
250                       digit->Strip(), digit->Counts(), counts, edep, mult));
251     
252     // Create a `RecPoint' on the output branch. 
253     AliFMDRecPoint* m = 
254       new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(), 
255                                             digit->Ring(), 
256                                             digit->Sector(),
257                                             digit->Strip(),
258                                             eta, phi, 
259                                             edep, mult);
260     (void)m; // Suppress warnings about unused variables. 
261     fNMult++;
262
263     fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(), 
264                              digit->Sector(),  digit->Strip(), mult);
265     fESDObj->SetEta(digit->Detector(), digit->Ring(), 
266                     digit->Sector(),  digit->Strip(), eta);
267   }
268 }
269
270 //____________________________________________________________________
271 UShort_t
272 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
273 {
274   // Member function to subtract the pedestal from a digit
275   // This implementation does nothing, but a derived class could over
276   // load this to subtract a pedestal that was given in a database or
277   // something like that. 
278
279   Int_t             counts = 0;
280   AliFMDParameters* param  = AliFMDParameters::Instance();
281   Float_t           pedM   = param->GetPedestal(digit->Detector(), 
282                                                 digit->Ring(), 
283                                                 digit->Sector(), 
284                                                 digit->Strip());
285   AliDebug(10, Form("Subtracting pedestal %f from signal %d", 
286                    pedM, digit->Counts()));
287   if (digit->Count3() > 0)      counts = digit->Count3();
288   else if (digit->Count2() > 0) counts = digit->Count2();
289   else                          counts = digit->Count1();
290   counts = TMath::Max(Int_t(counts - pedM), 0);
291   if (counts > 0) AliDebug(10, "Got a hit strip");
292   
293   return  UShort_t(counts);
294 }
295
296 //____________________________________________________________________
297 Float_t
298 AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit, 
299                                 Float_t      /* eta */, 
300                                 UShort_t     count) const
301 {
302   // Converts number of ADC counts to energy deposited. 
303   // Note, that this member function can be overloaded by derived
304   // classes to do strip-specific look-ups in databases or the like,
305   // to find the proper gain for a strip. 
306   // 
307   // In this simple version, we calculate the energy deposited as 
308   // 
309   //    EnergyDeposited = cos(theta) * gain * count
310   // 
311   // where 
312   // 
313   //           Pre_amp_MIP_Range
314   //    gain = ----------------- * Energy_deposited_per_MIP
315   //           ADC_channel_size    
316   // 
317   // is constant and the same for all strips.
318
319   // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
320   // Double_t edep  = TMath::Abs(TMath::Cos(theta)) * fGain * count;
321   AliFMDParameters* param = AliFMDParameters::Instance();
322   Float_t           gain  = param->GetPulseGain(digit->Detector(), 
323                                                 digit->Ring(), 
324                                                 digit->Sector(), 
325                                                 digit->Strip());
326   Double_t          edep  = count * gain;
327   AliDebug(10, Form("Converting counts %d to energy via factor %f", 
328                     count, gain));
329   return edep;
330 }
331
332 //____________________________________________________________________
333 Float_t
334 AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */, 
335                                          Float_t      edep) const
336 {
337   // Converts an energy signal to number of particles. 
338   // Note, that this member function can be overloaded by derived
339   // classes to do strip-specific look-ups in databases or the like,
340   // to find the proper gain for a strip. 
341   // 
342   // In this simple version, we calculate the multiplicity as 
343   // 
344   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
345   // 
346   // where 
347   //
348   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
349   // 
350   // is constant and the same for all strips 
351   AliFMDParameters* param   = AliFMDParameters::Instance();
352   Double_t          edepMIP = param->GetEdepMip();
353   Float_t           mult    = edep / edepMIP;
354   if (edep > 0) 
355     AliDebug(10, Form("Translating energy %f to multiplicity via "
356                      "divider %f->%f", edep, edepMIP, mult));
357   return mult;
358 }
359
360 //____________________________________________________________________
361 void
362 AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit, 
363                                          Float_t& eta, 
364                                          Float_t& phi) const
365 {
366   // Get the eta and phi of a digit 
367   // 
368   // Get geometry. 
369   AliFMDGeometry* geom = AliFMDGeometry::Instance();
370   Double_t x, y, z, r, theta;
371   geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(), 
372                     digit->Strip(), x, y, z);
373   // Correct for vertex offset. 
374   z     += fCurrentVertex;
375   phi   =  TMath::ATan2(y, x);
376   r     =  TMath::Sqrt(y * y + x * x);
377   theta =  TMath::ATan2(r, z);
378   eta   = -TMath::Log(TMath::Tan(theta / 2));
379 }
380
381       
382
383 //____________________________________________________________________
384 void 
385 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
386                              TTree*  /* clusterTree */,
387                              AliESD* esd) const
388 {
389   // nothing to be done
390   // FIXME: The vertex may not be known when Reconstruct is executed,
391   // so we may have to move some of that member function here. 
392   AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
393   // fESDObj->Print();
394
395   if (esd) { 
396     AliDebug(1, Form("Writing FMD data to ESD tree"));
397     esd->SetFMDData(fESDObj);
398   }
399 }
400
401
402 //____________________________________________________________________
403 void 
404 AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const 
405 {
406   // Cannot be used.  See member function with same name but with 2
407   // TTree arguments.   Make sure you do local reconstrucion 
408   AliDebug(1, Form("Calling FillESD with loader and tree"));
409   AliError("MayNotUse");
410 }
411 //____________________________________________________________________
412 void 
413 AliFMDReconstructor::Reconstruct(AliRunLoader*) const 
414 {
415   // Cannot be used.  See member function with same name but with 2
416   // TTree arguments.   Make sure you do local reconstrucion 
417   AliDebug(1, Form("Calling FillESD with loader"));
418   AliError("MayNotUse");
419 }
420 //____________________________________________________________________
421 void 
422 AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const 
423 {
424   // Cannot be used.  See member function with same name but with 2
425   // TTree arguments.   Make sure you do local reconstrucion 
426   AliDebug(1, Form("Calling FillESD with loader and raw reader"));
427   AliError("MayNotUse");
428 }
429 //____________________________________________________________________
430 void 
431 AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const 
432 {
433   // Cannot be used.  See member function with same name but with 2
434   // TTree arguments.   Make sure you do local reconstrucion 
435   AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
436   AliError("MayNotUse");
437 }
438 //____________________________________________________________________
439 void 
440 AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
441 {
442   // Cannot be used.  See member function with same name but with 2
443   // TTree arguments.   Make sure you do local reconstrucion 
444   AliDebug(1, Form("Calling FillESD with loader and ESD"));
445   AliError("MayNotUse");
446 }
447 //____________________________________________________________________
448 void 
449 AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const 
450 {
451   // Cannot be used.  See member function with same name but with 2
452   // TTree arguments.   Make sure you do local reconstrucion 
453   AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
454   AliError("MayNotUse");
455 }
456
457 //____________________________________________________________________
458 //
459 // EOF
460 //