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