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