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