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