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