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