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