Removed generated files@
[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 "AliFMDDebug.h"
38 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
39 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
40 #include "AliFMDAltroMapping.h"            // ALIFMDALTROMAPPING_H
41 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
42 #include "AliFMDSDigit.h"                  // ALIFMDDIGIT_H
43 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
44 #include "AliFMDRecoParam.h"               // ALIFMDRECOPARAM_H
45 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
46 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
47 #include "AliESDEvent.h"                   // ALIESDEVENT_H
48 #include "AliESDVertex.h"                  // ALIESDVERTEX_H
49 #include "AliESDTZERO.h"                   // ALIESDVERTEX_H
50 #include <AliESDFMD.h>                     // ALIESDFMD_H
51 #include <TMath.h>
52 #include <TH1.h>
53 #include <TH2.h>
54 #include <TFile.h>
55 #include <climits>
56 #include "AliFMDESDRevertexer.h"
57
58
59 class AliRawReader;
60
61 //____________________________________________________________________
62 ClassImp(AliFMDReconstructor)
63 #if 0
64   ; // This is here to keep Emacs for indenting the next line
65 #endif
66
67 //____________________________________________________________________
68 AliFMDReconstructor::AliFMDReconstructor() 
69   : AliReconstructor(),
70     fMult(0x0), 
71     fNMult(0),
72     fTreeR(0x0),
73     fCurrentVertex(0),
74     fESDObj(0x0),
75     fNoiseFactor(0),
76     fAngleCorrect(kTRUE),
77     fVertexType(kNoVertex),
78     fESD(0x0),
79     fDiagnostics(kFALSE),
80     fDiagStep1(0), 
81     fDiagStep2(0),
82     fDiagStep3(0),
83     fDiagStep4(0),
84     fDiagAll(0)
85 {
86   // Make a new FMD reconstructor object - default CTOR.  
87   SetNoiseFactor();
88   SetAngleCorrect();
89   if (AliDebugLevel() > 0) fDiagnostics = kTRUE;
90   for(Int_t det = 1; det<=3; det++) {
91     fZS[det-1]       = kFALSE;
92     fZSFactor[det-1] = 0;
93   }
94 }
95
96 //____________________________________________________________________
97 AliFMDReconstructor::~AliFMDReconstructor() 
98 {
99   // Destructor 
100   if (fMult)   fMult->Delete();
101   if (fMult)   delete fMult;
102   if (fESDObj) delete fESDObj;
103 }
104
105 //____________________________________________________________________
106 void 
107 AliFMDReconstructor::Init() 
108 {
109   // Initialize the reconstructor 
110
111   // Initialize the geometry 
112   AliFMDGeometry* geom = AliFMDGeometry::Instance();
113   geom->Init();
114   geom->InitTransformations();
115
116   // Initialize the parameters
117   AliFMDParameters* param = AliFMDParameters::Instance();
118   param->Init();
119   
120   // Current vertex position
121   fCurrentVertex = 0;
122   // Create array of reconstructed strip multiplicities 
123   // fMult = new TClonesArray("AliFMDRecPoint", 51200);
124   // Create ESD output object 
125   fESDObj = new AliESDFMD;
126   
127   // Check if we need diagnostics histograms 
128   if (!fDiagnostics) return;
129   AliInfo("Making diagnostics histograms");
130   fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
131                         1024, -.5, 1023.5, 1024, -.5, 1023.5);
132   fDiagStep1->SetDirectory(0);
133   fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
134   fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)", 
135                                         fNoiseFactor));
136   fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
137                         1024, -.5, 1023.5, 100, 0, 2);
138   fDiagStep2->SetDirectory(0);
139   fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
140   fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
141   fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
142                         100, 0., 2., 100, 0., 2.);
143   fDiagStep3->SetDirectory(0);
144   fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
145   fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
146   fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced", 
147                         100, 0., 2., 100, -.1, 19.9);
148   fDiagStep4->SetDirectory(0);
149   fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
150   fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
151   fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced", 
152                          1024, -.5, 1023.5, 100, -.1, 19.9);
153   fDiagAll->SetDirectory(0);
154   fDiagAll->GetXaxis()->SetTitle("ADC (read)");
155   fDiagAll->GetYaxis()->SetTitle("Multiplicity");
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   AliFMDDebug(1, ("Reading raw data into digits tree"));
165   if (!digitsTree) { 
166     AliError("No digits tree passed");
167     return;
168   }
169   static TClonesArray* array = new TClonesArray("AliFMDDigit");
170   digitsTree->Branch("FMD", &array);
171   array->Clear();
172   
173   AliFMDRawReader rawRead(reader, digitsTree);
174   // rawRead.SetSampleRate(fFMD->GetSampleRate());
175   // rawRead.Exec();
176   rawRead.ReadAdcs(array);
177
178   Int_t nWrite = digitsTree->Fill();
179   AliFMDDebug(1, ("Got a grand total of %d digits, wrote %d bytes to tree", 
180                    array->GetEntriesFast(), nWrite));
181
182   
183   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
184   for (size_t i = 1; i <= 3; i++) { 
185     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
186     fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
187   }
188 }
189
190 //____________________________________________________________________
191 void 
192 AliFMDReconstructor::GetVertex(AliESDEvent* esd) const
193 {
194   // Return the vertex to use. 
195   // This is obtained from the ESD object. 
196   // If not found, a warning is issued.
197   fVertexType    = kNoVertex;
198   fCurrentVertex = 0;
199   if (!esd) return;
200   
201   const AliESDVertex* vertex = esd->GetPrimaryVertex();
202   if (!vertex)        vertex = esd->GetPrimaryVertexSPD();
203   if (!vertex)        vertex = esd->GetPrimaryVertexTPC();
204   if (!vertex)        vertex = esd->GetVertex();
205
206   if (vertex) {
207     AliFMDDebug(2, ("Got %s (%s) from ESD: %f", 
208                     vertex->GetName(), vertex->GetTitle(), vertex->GetZv()));
209     fCurrentVertex = vertex->GetZv();
210     fVertexType    = kESDVertex;
211     return;
212   }
213   else if (esd->GetESDTZERO()) { 
214     AliFMDDebug(2, ("Got primary vertex from T0: %f", esd->GetT0zVertex()));
215     fCurrentVertex = esd->GetT0zVertex();
216     fVertexType    = kESDVertex;
217     return;
218   }
219   AliWarning("Didn't get any vertex from ESD or generator");
220 }
221   
222 //____________________________________________________________________
223 Int_t
224 AliFMDReconstructor::GetIdentifier() const
225 {
226   // Get the detector identifier. 
227   // Note the actual value is cached so that we do not 
228   // need to do many expensive string comparisons. 
229   static Int_t idx = AliReconstruction::GetDetIndex(GetDetectorName());
230   return idx;
231 }
232
233 //____________________________________________________________________
234 const AliFMDRecoParam*
235 AliFMDReconstructor::GetParameters() const
236 {
237   // Get the reconstruction parameters. 
238   // 
239   // Return: 
240   //   Pointer to reconstruction parameters or null if not found or wrong type
241   Int_t iDet = GetIdentifier(); // Was 12 - but changed on Cvetans request
242   const AliDetectorRecoParam* params = AliReconstructor::GetRecoParam(iDet);
243   if (!params || params->IsA() != AliFMDRecoParam::Class()) return 0;
244   return static_cast<const AliFMDRecoParam*>(params);
245 }
246
247 //____________________________________________________________________
248 void 
249 AliFMDReconstructor::UseRecoParam(Bool_t set) const
250 {
251   static Float_t savedNoiseFactor  = fNoiseFactor;
252   static Bool_t  savedAngleCorrect = fAngleCorrect;
253   if (set) { 
254     const AliFMDRecoParam* params  = GetParameters();
255     if (params) { 
256       fNoiseFactor  = params->NoiseFactor();
257       fAngleCorrect = params->AngleCorrect();
258     }
259     return;
260   }
261   fNoiseFactor  = savedNoiseFactor;
262   fAngleCorrect = savedAngleCorrect;
263 }
264   
265   
266 //____________________________________________________________________
267 void 
268 AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
269 {
270   // Loop over all entries of the ESD and mark 
271   // those that are dead as such 
272   // - otherwise put in the zero signal. 
273   AliFMDParameters* param = AliFMDParameters::Instance();
274
275   for (UShort_t d = 1; d <= 3; d++) { 
276     UShort_t nR = (d == 1 ? 1 : 2);
277
278     for (UShort_t q = 0; q < nR; q++) {
279       Char_t   r  = (q == 0 ? 'I' : 'O');
280       UShort_t nS = (q == 0 ?  20 :  40);
281       UShort_t nT = (q == 0 ? 512 : 256);
282
283       for (UShort_t s = 0; s < nS; s++) { 
284         for (UShort_t t = 0; t < nT; t++) {
285           if (param->IsDead(d, r, s, t)) { 
286             AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as dead", d, r, s, t));
287             esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
288             esd->SetEta(d, r, s, t, AliESDFMD::kInvalidEta);
289           }
290           else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) {
291             AliDebug(10, Form("Setting null signal in FMD%d%c[%2d,%3d]", d, r, s, t));
292             esd->SetMultiplicity(d, r, s, t, 0);
293           }
294         }
295       }
296     }
297   }
298 }
299
300 //____________________________________________________________________
301 void 
302 AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
303 {
304   // Reconstruct directly from raw data (no intermediate output on
305   // digit tree or rec point tree).  
306   // 
307   // Parameters: 
308   //   reader   Raw event reader 
309   //   ctree    Not used - 'cluster tree' to store rec-points in. 
310   AliFMDDebug(1, ("Reconstructing from raw reader"));
311   AliFMDRawReader rawReader(reader, 0);
312
313   UShort_t det, sec, str, fac;
314   Short_t  adc, oldDet = -1;
315   Bool_t   zs;
316   Char_t   rng;
317  
318   UseRecoParam(kTRUE);
319   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) { 
320     if (det != oldDet) { 
321       fZS[det-1]       = zs;
322       fZSFactor[det-1] = fac;
323       oldDet           = det;
324     }
325     ProcessSignal(det, rng, sec, str, adc);
326   }
327   UseRecoParam(kFALSE);
328   
329 }
330
331 //____________________________________________________________________
332 void 
333 AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
334 {
335   // Reconstruct directly from raw data (no intermediate output on
336   // digit tree or rec point tree).  
337   // 
338   // Parameters: 
339   //   reader   Raw event reader 
340   //   ctree    Not used. 
341   AliFMDRawReader rawReader(reader, 0);
342
343   UShort_t det, sec, str, sam, rat, fac;
344   Short_t  adc, oldDet = -1;
345   Bool_t   zs;
346   Char_t   rng;
347     
348   UseRecoParam(kTRUE);
349   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) { 
350     if (!rawReader.SelectSample(sam, rat)) continue;
351     if (det != oldDet) { 
352       fZS[det-1]       = zs;
353       fZSFactor[det-1] = fac;
354       oldDet           = det;
355     }
356     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
357   }
358   UseRecoParam(kFALSE);
359 }
360
361 //____________________________________________________________________
362 void 
363 AliFMDReconstructor::Reconstruct(TTree* digitsTree, 
364                                  TTree* clusterTree) const 
365 {
366   // Reconstruct event from digits in tree 
367   // Get the FMD branch holding the digits. 
368   // FIXME: The vertex may not be known yet, so we may have to move
369   // some of this to FillESD. 
370   // 
371   // Parameters: 
372   //   digitsTree       Pointer to a tree containing digits 
373   //   clusterTree      Pointer to output tree 
374   // 
375   if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
376
377   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
378   GetVertex(fESD);
379
380
381
382   static TClonesArray* digits = new TClonesArray("AliFMDDigit");
383   TBranch*      digitBranch   = digitsTree->GetBranch("FMD");
384   if (!digitBranch) {
385     Error("Exec", "No digit branch for the FMD found");
386     return;
387   }
388   digitBranch->SetAddress(&digits);
389
390   if (digits)  digits->Clear();
391   if (fMult)   fMult->Clear();
392   if (fESDObj) fESDObj->Clear();
393   
394   fNMult = 0;
395   fTreeR = clusterTree;
396   fTreeR->Branch("FMD", &fMult);
397   
398   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
399   digitBranch->GetEntry(0);
400   
401   AliFMDDebug(5, ("Processing digits"));
402   UseRecoParam(kTRUE);
403   ProcessDigits(digits);
404   UseRecoParam(kFALSE);
405
406   Int_t written = clusterTree->Fill();
407   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
408   // digits->Delete();
409   // delete digits;
410 }
411  
412
413 //____________________________________________________________________
414 void
415 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
416 {
417   // For each digit, find the pseudo rapdity, azimuthal angle, and
418   // number of corrected ADC counts, and pass it on to the algorithms
419   // used. 
420   // 
421   // Parameters: 
422   //    digits  Array of digits
423   // 
424   Int_t nDigits = digits->GetEntries();
425   AliFMDDebug(2, ("Got %d digits", nDigits));
426   fESDObj->SetNoiseFactor(fNoiseFactor);
427   fESDObj->SetAngleCorrected(fAngleCorrect);
428   for (Int_t i = 0; i < nDigits; i++) {
429     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
430     if (!digit) continue;
431     ProcessDigit(digit);
432   }
433 }
434
435 //____________________________________________________________________
436 void
437 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
438 {
439   UShort_t det = digit->Detector();
440   Char_t   rng = digit->Ring();
441   UShort_t sec = digit->Sector();
442   UShort_t str = digit->Strip();
443   Short_t  adc = digit->Counts();
444  
445   ProcessSignal(det, rng, sec, str, adc);
446 }
447
448 //____________________________________________________________________
449 void
450 AliFMDReconstructor::ProcessSignal(UShort_t det, 
451                                    Char_t   rng, 
452                                    UShort_t sec, 
453                                    UShort_t str, 
454                                    Short_t  adc) const
455 {
456   // Process the signal from a single strip 
457   // 
458   // Parameters: 
459   //    det     Detector ID
460   //    rng     Ring ID
461   //    sec     Sector ID
462   //    rng     Strip ID
463   //    adc     ADC counts
464   // 
465   AliFMDParameters* param  = AliFMDParameters::Instance();
466   // Check that the strip is not marked as dead 
467   if (param->IsDead(det, rng, sec, str)) {
468     AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
469     return;
470   }
471   
472   // digit->Print();
473   // Get eta and phi 
474   Float_t eta, phi;
475   PhysicalCoordinates(det, rng, sec, str, eta, phi);
476     
477   // Substract pedestal. 
478   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
479   if(counts == USHRT_MAX) return;
480   
481     // Gain match digits. 
482   Double_t edep     = Adc2Energy(det, rng, sec, str, eta, counts);
483   // Get rid of nonsense energy
484   if(edep < 0)  return;
485   
486   // Make rough multiplicity 
487   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
488   // Get rid of nonsense mult
489   //if (mult > 20) { 
490   //  AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 "
491   //                "(ADC: %d, Energy: %f)", det, rng, sec, str, mult, 
492   //                counts, edep));
493   // }
494   if (mult < 0)  return; 
495   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
496                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
497                   det, rng, sec, str, adc, counts, edep, mult));
498   
499   // Create a `RecPoint' on the output branch. 
500   if (fMult) {
501     AliFMDRecPoint* m = 
502       new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str, 
503                                             eta, phi, edep, mult);
504     (void)m; // Suppress warnings about unused variables. 
505     fNMult++;
506   }
507   
508   fESDObj->SetMultiplicity(det, rng, sec, str, mult);
509   fESDObj->SetEta(det, rng, sec, str, eta);
510   
511   if (fDiagAll) fDiagAll->Fill(adc, mult);  
512
513 }
514
515 //____________________________________________________________________
516 void
517 AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits, 
518                                     UShort_t det, 
519                                     Char_t   rng, 
520                                     UShort_t sec, 
521                                     UShort_t str, 
522                                     UShort_t /* sam */,
523                                     Short_t  adc) const
524 {
525   // Process the signal from a single strip 
526   // 
527   // Parameters: 
528   //    det     Detector ID
529   //    rng     Ring ID
530   //    sec     Sector ID
531   //    rng     Strip ID
532   //    adc     ADC counts
533   // 
534   AliFMDParameters* param  = AliFMDParameters::Instance();
535   // Check that the strip is not marked as dead 
536   if (param->IsDead(det, rng, sec, str)) {
537     AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
538     return;
539   }
540   
541   // Substract pedestal. 
542   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
543   if(counts == USHRT_MAX || counts == 0) return;
544   
545     // Gain match digits. 
546   Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
547   // Get rid of nonsense energy
548   if(edep < 0)  return;
549
550   Int_t n = sdigits->GetEntriesFast();
551   // AliFMDSDigit* sdigit = 
552   new ((*sdigits)[n]) 
553     AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
554   // sdigit->SetCount(sam, counts);
555 }
556
557 //____________________________________________________________________
558 UShort_t
559 AliFMDReconstructor::SubtractPedestal(UShort_t det, 
560                                       Char_t   rng, 
561                                       UShort_t sec, 
562                                       UShort_t str, 
563                                       UShort_t adc, 
564                                       Float_t  noiseFactor,
565                                       Bool_t   zsEnabled, 
566                                       UShort_t zsNoiseFactor) const
567 {
568   AliFMDParameters* param  = AliFMDParameters::Instance();
569   Float_t           ped    = (zsEnabled ? 0 : 
570                                 param->GetPedestal(det, rng, sec, str));
571   Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
572   if(ped < 0 || noise < 0) { 
573     AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
574                          "for FMD%d%c[%02d,%03d]", 
575                     ped, noise, det, rng, sec, str));
576     return USHRT_MAX;
577   }
578   AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
579                          "(%s w/factor %d, noise factor %f, "
580                          "pedestal %8.2f+/-%8.2f)",
581                          det, rng, sec, str, adc, 
582                          (zsEnabled ? "zs'ed" : "straight"), 
583                          zsNoiseFactor, noiseFactor, ped, noise));
584
585   Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
586   counts =  TMath::Max(Int_t(counts), 0);
587   // Calculate the noise factor for suppressing remenants of the noise
588   // peak.  If we have done on-line zero suppression, we only check
589   // for noise signals that are larger than the suppressed noise.  If
590   // the noise factor used on line is larger than the factor used
591   // here, we do not do this check at all.  
592   // 
593   // For example:
594   //    Online factor  |  Read factor |  Result 
595   //    ---------------+--------------+-------------------------------
596   //           2       |      3       | Check if signal > 1 * noise
597   //           3       |      3       | Check if signal > 0
598   //           3       |      2       | Check if signal > 0
599   //
600   // In this way, we make sure that we do not suppress away too much
601   // data, and that the read-factor is the most stringent cut. 
602   Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
603   if (counts < noise * nf) counts = 0;
604   if (counts > 0) AliDebugClass(15, "Got a hit strip");
605
606   UShort_t ret = counts < 0 ? 0 : counts;
607   return ret;
608 }
609
610
611 //____________________________________________________________________
612 UShort_t
613 AliFMDReconstructor::SubtractPedestal(UShort_t det, 
614                                       Char_t   rng, 
615                                       UShort_t sec, 
616                                       UShort_t str, 
617                                       Short_t  adc) const
618 {
619   // Member function to subtract the pedestal from a digit
620   //
621   // Parameters: 
622   //    det     Detector ID
623   //    rng     Ring ID
624   //    sec     Sector ID
625   //    rng     Strip ID
626   //    adc     # of ADC counts
627   // Return:
628   //    Pedestal subtracted signal or USHRT_MAX in case of problems 
629   //
630   UShort_t counts = SubtractPedestal(det, rng, sec, str, adc, 
631                                      fNoiseFactor, fZS[det-1], 
632                                      fZSFactor[det-1]);
633   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
634   
635   return counts;
636 }
637
638 //____________________________________________________________________
639 Float_t
640 AliFMDReconstructor::Adc2Energy(UShort_t det, 
641                                 Char_t   rng, 
642                                 UShort_t sec, 
643                                 UShort_t str, 
644                                 UShort_t count) const
645 {
646   // Converts number of ADC counts to energy deposited. 
647   // Note, that this member function can be overloaded by derived
648   // classes to do strip-specific look-ups in databases or the like,
649   // to find the proper gain for a strip. 
650   // 
651   // In the first simple version, we calculate the energy deposited as 
652   // 
653   //    EnergyDeposited = cos(theta) * gain * count
654   // 
655   // where 
656   // 
657   //           Pre_amp_MIP_Range
658   //    gain = ----------------- * Energy_deposited_per_MIP
659   //           ADC_channel_size    
660   // 
661   // is constant and the same for all strips.
662   //
663   // For the production we use the conversion measured in the NBI lab.
664   // The total conversion is then:
665   // 
666   //    gain = ADC / DAC
667   // 
668   //                  EdepMip * count
669   //      => energy = ----------------
670   //                  gain * DACPerADC
671   // 
672   // Parameters: 
673   //    det     Detector ID
674   //    rng     Ring ID
675   //    sec     Sector ID
676   //    rng     Strip ID
677   //    counts  Number of ADC counts over pedestal
678   // Return 
679   //    The energy deposited in a single strip, or -1 in case of problems
680   //
681   if (count <= 0) return 0;
682   AliFMDParameters* param = AliFMDParameters::Instance();
683   Float_t           gain  = param->GetPulseGain(det, rng, sec, str);
684   // 'Tagging' bad gains as bad energy
685   if (gain < 0) { 
686     AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]", 
687                     gain, det, rng, sec, str));
688     return -1;
689   }
690   AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)", 
691                   count, gain,param->GetDACPerMIP()));
692
693   Double_t edep  = ((count * param->GetEdepMip()) 
694                     / (gain * param->GetDACPerMIP()));
695   return edep;
696 }
697
698 //____________________________________________________________________
699 Float_t
700 AliFMDReconstructor::Adc2Energy(UShort_t det, 
701                                 Char_t   rng, 
702                                 UShort_t sec, 
703                                 UShort_t str, 
704                                 Float_t  eta, 
705                                 UShort_t count) const
706 {
707   // Converts number of ADC counts to energy deposited. 
708   // Note, that this member function can be overloaded by derived
709   // classes to do strip-specific look-ups in databases or the like,
710   // to find the proper gain for a strip. 
711   // 
712   // In the first simple version, we calculate the energy deposited as 
713   // 
714   //    EnergyDeposited = cos(theta) * gain * count
715   // 
716   // where 
717   // 
718   //           Pre_amp_MIP_Range
719   //    gain = ----------------- * Energy_deposited_per_MIP
720   //           ADC_channel_size    
721   // 
722   // is constant and the same for all strips.
723   //
724   // For the production we use the conversion measured in the NBI lab.
725   // The total conversion is then:
726   // 
727   //    gain = ADC / DAC
728   // 
729   //                  EdepMip * count
730   //      => energy = ----------------
731   //                  gain * DACPerADC
732   // 
733   // Parameters: 
734   //    det     Detector ID
735   //    rng     Ring ID
736   //    sec     Sector ID
737   //    rng     Strip ID
738   //    eta     Psuedo-rapidity
739   //    counts  Number of ADC counts over pedestal
740   // Return 
741   //    The energy deposited in a single strip, or -1 in case of problems
742   //
743   Double_t edep = Adc2Energy(det, rng, sec, str, count);
744   
745   if (fDiagStep2) fDiagStep2->Fill(count, edep);  
746   if (fAngleCorrect) {
747     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
748     Double_t corr  = TMath::Abs(TMath::Cos(theta));
749     Double_t cedep = corr * edep;
750     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
751                       edep, corr, cedep, eta, theta));
752     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);  
753     edep = cedep;
754   }
755   return edep;
756 }
757
758 //____________________________________________________________________
759 Float_t
760 AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/, 
761                                          Char_t   /*rng*/, 
762                                          UShort_t /*sec*/, 
763                                          UShort_t /*str*/, 
764                                          Float_t  edep) const
765 {
766   // Converts an energy signal to number of particles. 
767   // Note, that this member function can be overloaded by derived
768   // classes to do strip-specific look-ups in databases or the like,
769   // to find the proper gain for a strip. 
770   // 
771   // In this simple version, we calculate the multiplicity as 
772   // 
773   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
774   // 
775   // where 
776   //
777   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
778   // 
779   // is constant and the same for all strips 
780   //
781   // Parameters: 
782   //    det     Detector ID
783   //    rng     Ring ID
784   //    sec     Sector ID
785   //    rng     Strip ID
786   //    edep    Energy deposited in a single strip
787   // Return 
788   //    The "bare" multiplicity corresponding to the energy deposited
789   AliFMDParameters* param   = AliFMDParameters::Instance();
790   Double_t          edepMIP = param->GetEdepMip();
791   Float_t           mult    = edep / edepMIP;
792 #if 0
793   if (edep > 0) 
794     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
795                      "divider %f->%f", edep, edepMIP, mult));
796 #endif
797   if (fDiagStep4) fDiagStep4->Fill(edep, mult);  
798   return mult;
799 }
800
801 //____________________________________________________________________
802 void
803 AliFMDReconstructor::PhysicalCoordinates(UShort_t det, 
804                                          Char_t   rng, 
805                                          UShort_t sec, 
806                                          UShort_t str, 
807                                          Float_t& eta, 
808                                          Float_t& phi) const
809 {
810   // Get the eta and phi of a digit 
811   // 
812   // Parameters: 
813   //    det     Detector ID
814   //    rng     Ring ID
815   //    sec     Sector ID
816   //    rng     Strip ID
817   //    eta     On return, contains the psuedo-rapidity of the strip
818   //    phi     On return, contains the azimuthal angle of the strip
819   // 
820   AliFMDGeometry* geom = AliFMDGeometry::Instance();
821   Double_t x, y, z, r, theta, deta, dphi;
822   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
823
824   // Correct for vertex offset. 
825   z     -= fCurrentVertex;
826   AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
827   eta = deta;
828   phi = dphi;
829 }
830
831 namespace { 
832   class ESDPrinter : public AliESDFMD::ForOne
833   {
834   public:
835     ESDPrinter() {}
836     Bool_t operator()(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
837                       Float_t m, Float_t e)
838     {
839       if (m > 0 && m != AliESDFMD::kInvalidMult) 
840         printf("  FMD%d%c[%2d,%3d] = %6.3f / %6.3f\n", d, r, s, t, m, e);
841       return kTRUE;
842     }
843   };
844 }
845
846
847 //____________________________________________________________________
848 void 
849 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
850                              TTree*  /* clusterTree */,
851                              AliESDEvent* esd) const
852 {
853   // nothing to be done
854   // FIXME: The vertex may not be known when Reconstruct is executed,
855   // so we may have to move some of that member function here. 
856   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
857   // fESDObj->Print();
858
859   // Fix up ESD so that only truely dead channels get the kInvalidMult flag. 
860   MarkDeadChannels(fESDObj);
861
862   Double_t oldVz = fCurrentVertex;
863   GetVertex(esd);
864   if (fVertexType != kNoVertex) { 
865     AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
866                     fCurrentVertex, oldVz));
867     AliFMDESDRevertexer revertexer;
868     revertexer.Revertex(fESDObj, fCurrentVertex);
869   }
870   
871   if (AliDebugLevel() > 10) { 
872     ESDPrinter p;
873     fESDObj->ForEach(p);
874   }
875
876   if (esd) { 
877     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
878     esd->SetFMDData(fESDObj);
879   }
880
881   if (!fDiagnostics || !esd) return;
882   static bool first = true;
883   // This is most likely NOT the event number you'd like to use. It
884   // has nothing to do with the 'real' event number. 
885   // - That's OK.  We just use it for the name of the directory -
886   // nothing else.  Christian
887   Int_t evno = esd->GetEventNumberInFile(); 
888   AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
889   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
890   first = false;
891   f.cd(); 
892   TDirectory* d = f.mkdir(Form("%03d", evno), 
893                           Form("Diagnostics histograms for event # %d", evno));
894   d->cd();
895   if (fDiagStep1) fDiagStep1->Write();
896   if (fDiagStep2) fDiagStep2->Write();
897   if (fDiagStep3) fDiagStep3->Write();
898   if (fDiagStep4) fDiagStep4->Write();
899   if (fDiagAll)   fDiagAll->Write();
900   d->Write();
901   f.Write();
902   f.Close();
903
904   if (fDiagStep1) fDiagStep1->Reset();
905   if (fDiagStep2) fDiagStep2->Reset();
906   if (fDiagStep3) fDiagStep3->Reset();
907   if (fDiagStep4) fDiagStep4->Reset();
908   if (fDiagAll)   fDiagAll->Reset();
909 }
910
911 //____________________________________________________________________
912 void
913 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree, 
914                              AliESDEvent* esd) const
915 {
916   TTree* dummy = 0;
917   FillESD(dummy, clusterTree, esd);
918 }
919 //____________________________________________________________________
920 //
921 // EOF
922 //