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