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