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