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