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