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