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