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