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