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