Added the offline trigger functions for use with the first physics.
[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
ddaa8027 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.
3991ca94 304 AliFMDDebug(1, ("Reconstructing from raw reader"));
50b9d194 305 AliFMDRawReader rawReader(reader, 0);
306
307 UShort_t det, sec, str, fac;
308 Short_t adc, oldDet = -1;
309 Bool_t zs;
310 Char_t rng;
818fff8d 311
312 UseRecoParam(kTRUE);
50b9d194 313 while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
314 if (det != oldDet) {
315 fZS[det-1] = zs;
316 fZSFactor[det-1] = fac;
317 oldDet = det;
318 }
319 ProcessSignal(det, rng, sec, str, adc);
320 }
818fff8d 321 UseRecoParam(kFALSE);
322
ddaa8027 323}
324
faf80567 325//____________________________________________________________________
326void
327AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
328{
329 // Reconstruct directly from raw data (no intermediate output on
330 // digit tree or rec point tree).
331 //
332 // Parameters:
333 // reader Raw event reader
334 // ctree Not used.
335 AliFMDRawReader rawReader(reader, 0);
336
337 UShort_t det, sec, str, sam, rat, fac;
338 Short_t adc, oldDet = -1;
339 Bool_t zs;
340 Char_t rng;
341
818fff8d 342 UseRecoParam(kTRUE);
faf80567 343 while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) {
344 if (!rawReader.SelectSample(sam, rat)) continue;
345 if (det != oldDet) {
346 fZS[det-1] = zs;
347 fZSFactor[det-1] = fac;
348 oldDet = det;
349 }
350 DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
351 }
818fff8d 352 UseRecoParam(kFALSE);
faf80567 353}
354
9684be2f 355//____________________________________________________________________
356void
357AliFMDReconstructor::Reconstruct(TTree* digitsTree,
358 TTree* clusterTree) const
359{
360 // Reconstruct event from digits in tree
361 // Get the FMD branch holding the digits.
362 // FIXME: The vertex may not be known yet, so we may have to move
363 // some of this to FillESD.
50b9d194 364 //
365 // Parameters:
366 // digitsTree Pointer to a tree containing digits
367 // clusterTree Pointer to output tree
368 //
3991ca94 369 AliFMDDebug(1, ("Reconstructing from digits in a tree"));
8983e5ae 370 GetVertex(fESD);
b04d1f8e 371
1a1fdef7 372 TBranch *digitBranch = digitsTree->GetBranch("FMD");
4347b38f 373 if (!digitBranch) {
1a1fdef7 374 Error("Exec", "No digit branch for the FMD found");
375 return;
4347b38f 376 }
0abe7182 377 TClonesArray* digits = new TClonesArray("AliFMDDigit");
1a1fdef7 378 digitBranch->SetAddress(&digits);
4347b38f 379
8f6ee336 380 if (fMult) fMult->Clear();
381 if (fESDObj) fESDObj->Clear();
382
383 fNMult = 0;
384 fTreeR = clusterTree;
385 fTreeR->Branch("FMD", &fMult);
386
f95a63c4 387 AliFMDDebug(5, ("Getting entry 0 from digit branch"));
1a1fdef7 388 digitBranch->GetEntry(0);
389
3991ca94 390 AliFMDDebug(5, ("Processing digits"));
818fff8d 391 UseRecoParam(kTRUE);
e802be3e 392 ProcessDigits(digits);
818fff8d 393 UseRecoParam(kFALSE);
e802be3e 394
8f6ee336 395 Int_t written = clusterTree->Fill();
f95a63c4 396 AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
0abe7182 397 digits->Delete();
398 delete digits;
4347b38f 399}
1a1fdef7 400
8f6ee336 401
4347b38f 402//____________________________________________________________________
e802be3e 403void
404AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f 405{
69b696b9 406 // For each digit, find the pseudo rapdity, azimuthal angle, and
407 // number of corrected ADC counts, and pass it on to the algorithms
408 // used.
50b9d194 409 //
410 // Parameters:
411 // digits Array of digits
412 //
e802be3e 413 Int_t nDigits = digits->GetEntries();
3991ca94 414 AliFMDDebug(2, ("Got %d digits", nDigits));
a9579262 415 fESDObj->SetNoiseFactor(fNoiseFactor);
416 fESDObj->SetAngleCorrected(fAngleCorrect);
e802be3e 417 for (Int_t i = 0; i < nDigits; i++) {
418 AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
50b9d194 419 if (!digit) continue;
420 ProcessDigit(digit);
421 }
422}
8f6ee336 423
50b9d194 424//____________________________________________________________________
425void
426AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
427{
428 UShort_t det = digit->Detector();
429 Char_t rng = digit->Ring();
430 UShort_t sec = digit->Sector();
431 UShort_t str = digit->Strip();
432 Short_t adc = digit->Counts();
433
434 ProcessSignal(det, rng, sec, str, adc);
435}
436
437//____________________________________________________________________
438void
439AliFMDReconstructor::ProcessSignal(UShort_t det,
440 Char_t rng,
441 UShort_t sec,
442 UShort_t str,
443 Short_t adc) const
444{
445 // Process the signal from a single strip
446 //
447 // Parameters:
448 // det Detector ID
449 // rng Ring ID
450 // sec Sector ID
451 // rng Strip ID
452 // adc ADC counts
453 //
454 AliFMDParameters* param = AliFMDParameters::Instance();
455 // Check that the strip is not marked as dead
456 if (param->IsDead(det, rng, sec, str)) {
3991ca94 457 AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
50b9d194 458 return;
459 }
460
461 // digit->Print();
462 // Get eta and phi
463 Float_t eta, phi;
464 PhysicalCoordinates(det, rng, sec, str, eta, phi);
4347b38f 465
50b9d194 466 // Substract pedestal.
467 UShort_t counts = SubtractPedestal(det, rng, sec, str, adc);
468 if(counts == USHRT_MAX) return;
469
8f6ee336 470 // Gain match digits.
50b9d194 471 Double_t edep = Adc2Energy(det, rng, sec, str, eta, counts);
472 // Get rid of nonsense energy
473 if(edep < 0) return;
474
475 // Make rough multiplicity
476 Double_t mult = Energy2Multiplicity(det, rng, sec, str, edep);
477 // Get rid of nonsense mult
478 if (mult < 0) return;
3991ca94 479 AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
50b9d194 480 "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
481 det, rng, sec, str, adc, counts, edep, mult));
482
483 // Create a `RecPoint' on the output branch.
484 if (fMult) {
485 AliFMDRecPoint* m =
486 new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
487 eta, phi, edep, mult);
488 (void)m; // Suppress warnings about unused variables.
489 fNMult++;
4347b38f 490 }
50b9d194 491
492 fESDObj->SetMultiplicity(det, rng, sec, str, mult);
493 fESDObj->SetEta(det, rng, sec, str, eta);
494
495 if (fDiagAll) fDiagAll->Fill(adc, mult);
496
4347b38f 497}
8f6ee336 498
faf80567 499//____________________________________________________________________
500void
501AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits,
502 UShort_t det,
503 Char_t rng,
504 UShort_t sec,
505 UShort_t str,
506 UShort_t /* sam */,
507 Short_t adc) const
508{
509 // Process the signal from a single strip
510 //
511 // Parameters:
512 // det Detector ID
513 // rng Ring ID
514 // sec Sector ID
515 // rng Strip ID
516 // adc ADC counts
517 //
518 AliFMDParameters* param = AliFMDParameters::Instance();
519 // Check that the strip is not marked as dead
520 if (param->IsDead(det, rng, sec, str)) {
521 AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
522 return;
523 }
524
525 // Substract pedestal.
526 UShort_t counts = SubtractPedestal(det, rng, sec, str, adc);
527 if(counts == USHRT_MAX || counts == 0) return;
528
529 // Gain match digits.
530 Double_t edep = Adc2Energy(det, rng, sec, str, counts);
531 // Get rid of nonsense energy
532 if(edep < 0) return;
533
534 Int_t n = sdigits->GetEntriesFast();
535 // AliFMDSDigit* sdigit =
536 new ((*sdigits)[n])
537 AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
538 // sdigit->SetCount(sam, counts);
539}
540
541//____________________________________________________________________
542UShort_t
543AliFMDReconstructor::SubtractPedestal(UShort_t det,
544 Char_t rng,
545 UShort_t sec,
546 UShort_t str,
547 UShort_t adc,
548 Float_t noiseFactor,
549 Bool_t zsEnabled,
550 UShort_t zsNoiseFactor) const
551{
552 AliFMDParameters* param = AliFMDParameters::Instance();
553 Float_t ped = (zsEnabled ? 0 :
554 param->GetPedestal(det, rng, sec, str));
555 Float_t noise = param->GetPedestalWidth(det, rng, sec, str);
556 if(ped < 0 || noise < 0) {
557 AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
558 "for FMD%d%c[%02d,%03d]",
559 ped, noise, det, rng, sec, str));
560 return USHRT_MAX;
561 }
562 AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
563 "(%s w/factor %d, noise factor %f, "
564 "pedestal %8.2f+/-%8.2f)",
565 det, rng, sec, str, adc,
566 (zsEnabled ? "zs'ed" : "straight"),
567 zsNoiseFactor, noiseFactor, ped, noise));
568
569 Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
570 counts = TMath::Max(Int_t(counts), 0);
571 // Calculate the noise factor for suppressing remenants of the noise
572 // peak. If we have done on-line zero suppression, we only check
573 // for noise signals that are larger than the suppressed noise. If
574 // the noise factor used on line is larger than the factor used
575 // here, we do not do this check at all.
576 //
577 // For example:
578 // Online factor | Read factor | Result
579 // ---------------+--------------+-------------------------------
580 // 2 | 3 | Check if signal > 1 * noise
581 // 3 | 3 | Check if signal > 0
582 // 3 | 2 | Check if signal > 0
583 //
584 // In this way, we make sure that we do not suppress away too much
585 // data, and that the read-factor is the most stringent cut.
586 Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
587 if (counts < noise * nf) counts = 0;
588 if (counts > 0) AliDebugClass(15, "Got a hit strip");
589
590 return counts;
591}
50b9d194 592
593
1a1fdef7 594//____________________________________________________________________
595UShort_t
50b9d194 596AliFMDReconstructor::SubtractPedestal(UShort_t det,
597 Char_t rng,
598 UShort_t sec,
599 UShort_t str,
600 Short_t adc) const
1a1fdef7 601{
602 // Member function to subtract the pedestal from a digit
50b9d194 603 //
604 // Parameters:
605 // det Detector ID
606 // rng Ring ID
607 // sec Sector ID
608 // rng Strip ID
609 // adc # of ADC counts
610 // Return:
611 // Pedestal subtracted signal or USHRT_MAX in case of problems
612 //
faf80567 613 UShort_t counts = SubtractPedestal(det, rng, sec, str, adc,
614 fNoiseFactor, fZS[det-1],
615 fZSFactor[det-1]);
9684be2f 616 if (fDiagStep1) fDiagStep1->Fill(adc, counts);
8f6ee336 617
faf80567 618 return counts;
1a1fdef7 619}
620
8f6ee336 621//____________________________________________________________________
622Float_t
50b9d194 623AliFMDReconstructor::Adc2Energy(UShort_t det,
624 Char_t rng,
625 UShort_t sec,
626 UShort_t str,
50b9d194 627 UShort_t count) const
8f6ee336 628{
629 // Converts number of ADC counts to energy deposited.
630 // Note, that this member function can be overloaded by derived
631 // classes to do strip-specific look-ups in databases or the like,
632 // to find the proper gain for a strip.
633 //
68aba90a 634 // In the first simple version, we calculate the energy deposited as
8f6ee336 635 //
636 // EnergyDeposited = cos(theta) * gain * count
637 //
638 // where
639 //
640 // Pre_amp_MIP_Range
641 // gain = ----------------- * Energy_deposited_per_MIP
642 // ADC_channel_size
643 //
644 // is constant and the same for all strips.
68aba90a 645 //
646 // For the production we use the conversion measured in the NBI lab.
647 // The total conversion is then:
50b9d194 648 //
649 // gain = ADC / DAC
650 //
651 // EdepMip * count
652 // => energy = ----------------
653 // gain * DACPerADC
654 //
655 // Parameters:
656 // det Detector ID
657 // rng Ring ID
658 // sec Sector ID
659 // rng Strip ID
50b9d194 660 // counts Number of ADC counts over pedestal
661 // Return
662 // The energy deposited in a single strip, or -1 in case of problems
68aba90a 663 //
a9579262 664 if (count <= 0) return 0;
8f6ee336 665 AliFMDParameters* param = AliFMDParameters::Instance();
50b9d194 666 Float_t gain = param->GetPulseGain(det, rng, sec, str);
68aba90a 667 // 'Tagging' bad gains as bad energy
668 if (gain < 0) {
50b9d194 669 AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
670 gain, det, rng, sec, str));
68aba90a 671 return -1;
672 }
50b9d194 673 AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
68aba90a 674 count, gain,param->GetDACPerMIP()));
a9579262 675
50b9d194 676 Double_t edep = ((count * param->GetEdepMip())
677 / (gain * param->GetDACPerMIP()));
faf80567 678 return edep;
679}
680
681//____________________________________________________________________
682Float_t
683AliFMDReconstructor::Adc2Energy(UShort_t det,
684 Char_t rng,
685 UShort_t sec,
686 UShort_t str,
687 Float_t eta,
688 UShort_t count) const
689{
690 // Converts number of ADC counts to energy deposited.
691 // Note, that this member function can be overloaded by derived
692 // classes to do strip-specific look-ups in databases or the like,
693 // to find the proper gain for a strip.
694 //
695 // In the first simple version, we calculate the energy deposited as
696 //
697 // EnergyDeposited = cos(theta) * gain * count
698 //
699 // where
700 //
701 // Pre_amp_MIP_Range
702 // gain = ----------------- * Energy_deposited_per_MIP
703 // ADC_channel_size
704 //
705 // is constant and the same for all strips.
706 //
707 // For the production we use the conversion measured in the NBI lab.
708 // The total conversion is then:
709 //
710 // gain = ADC / DAC
711 //
712 // EdepMip * count
713 // => energy = ----------------
714 // gain * DACPerADC
715 //
716 // Parameters:
717 // det Detector ID
718 // rng Ring ID
719 // sec Sector ID
720 // rng Strip ID
721 // eta Psuedo-rapidity
722 // counts Number of ADC counts over pedestal
723 // Return
724 // The energy deposited in a single strip, or -1 in case of problems
725 //
726 Double_t edep = Adc2Energy(det, rng, sec, str, count);
727
9684be2f 728 if (fDiagStep2) fDiagStep2->Fill(count, edep);
a9579262 729 if (fAngleCorrect) {
9684be2f 730 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
731 Double_t corr = TMath::Abs(TMath::Cos(theta));
732 Double_t cedep = corr * edep;
f95a63c4 733 AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
9684be2f 734 edep, corr, cedep, eta, theta));
735 if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
736 edep = cedep;
a9579262 737 }
8f6ee336 738 return edep;
739}
740
741//____________________________________________________________________
742Float_t
50b9d194 743AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
744 Char_t /*rng*/,
745 UShort_t /*sec*/,
746 UShort_t /*str*/,
747 Float_t edep) const
8f6ee336 748{
749 // Converts an energy signal to number of particles.
750 // Note, that this member function can be overloaded by derived
751 // classes to do strip-specific look-ups in databases or the like,
752 // to find the proper gain for a strip.
753 //
754 // In this simple version, we calculate the multiplicity as
755 //
756 // multiplicity = Energy_deposited / Energy_deposited_per_MIP
757 //
758 // where
759 //
760 // Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
761 //
762 // is constant and the same for all strips
50b9d194 763 //
764 // Parameters:
765 // det Detector ID
766 // rng Ring ID
767 // sec Sector ID
768 // rng Strip ID
769 // edep Energy deposited in a single strip
770 // Return
771 // The "bare" multiplicity corresponding to the energy deposited
8f6ee336 772 AliFMDParameters* param = AliFMDParameters::Instance();
773 Double_t edepMIP = param->GetEdepMip();
774 Float_t mult = edep / edepMIP;
68aba90a 775#if 0
8f6ee336 776 if (edep > 0)
f95a63c4 777 AliFMDDebug(15, ("Translating energy %f to multiplicity via "
8f6ee336 778 "divider %f->%f", edep, edepMIP, mult));
68aba90a 779#endif
9684be2f 780 if (fDiagStep4) fDiagStep4->Fill(edep, mult);
8f6ee336 781 return mult;
782}
783
784//____________________________________________________________________
785void
50b9d194 786AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
787 Char_t rng,
788 UShort_t sec,
789 UShort_t str,
8f6ee336 790 Float_t& eta,
791 Float_t& phi) const
792{
793 // Get the eta and phi of a digit
794 //
50b9d194 795 // Parameters:
796 // det Detector ID
797 // rng Ring ID
798 // sec Sector ID
799 // rng Strip ID
800 // eta On return, contains the psuedo-rapidity of the strip
801 // phi On return, contains the azimuthal angle of the strip
802 //
9b48326f 803 AliFMDGeometry* geom = AliFMDGeometry::Instance();
bf000c32 804 Double_t x, y, z, r, theta;
50b9d194 805 geom->Detector2XYZ(det, rng, sec, str, x, y, z);
bf000c32 806 // Correct for vertex offset.
8161c5f7 807 z -= fCurrentVertex;
bf000c32 808 phi = TMath::ATan2(y, x);
809 r = TMath::Sqrt(y * y + x * x);
810 theta = TMath::ATan2(r, z);
811 eta = -TMath::Log(TMath::Tan(theta / 2));
8f6ee336 812}
813
814
815
4347b38f 816//____________________________________________________________________
817void
1a1fdef7 818AliFMDReconstructor::FillESD(TTree* /* digitsTree */,
819 TTree* /* clusterTree */,
af885e0f 820 AliESDEvent* esd) const
121a60bd 821{
42403906 822 // nothing to be done
69b696b9 823 // FIXME: The vertex may not be known when Reconstruct is executed,
824 // so we may have to move some of that member function here.
f95a63c4 825 AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
8f6ee336 826 // fESDObj->Print();
827
50b9d194 828 Double_t oldVz = fCurrentVertex;
8983e5ae 829 GetVertex(esd);
50b9d194 830 if (fVertexType != kNoVertex) {
831 AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
832 fCurrentVertex, oldVz));
833 AliFMDESDRevertexer revertexer;
834 revertexer.Revertex(fESDObj, fCurrentVertex);
835 }
836
8f6ee336 837 if (esd) {
f95a63c4 838 AliFMDDebug(2, ("Writing FMD data to ESD tree"));
8f6ee336 839 esd->SetFMDData(fESDObj);
a3537838 840 }
9684be2f 841
842 if (!fDiagnostics || !esd) return;
843 static bool first = true;
69893a66 844 // This is most likely NOT the event number you'd like to use. It
845 // has nothing to do with the 'real' event number.
846 // - That's OK. We just use it for the name of the directory -
847 // nothing else. Christian
848 Int_t evno = esd->GetEventNumberInFile();
3991ca94 849 AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
9684be2f 850 TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
851 first = false;
852 f.cd();
853 TDirectory* d = f.mkdir(Form("%03d", evno),
854 Form("Diagnostics histograms for event # %d", evno));
855 d->cd();
856 if (fDiagStep1) fDiagStep1->Write();
857 if (fDiagStep2) fDiagStep2->Write();
858 if (fDiagStep3) fDiagStep3->Write();
859 if (fDiagStep4) fDiagStep4->Write();
860 if (fDiagAll) fDiagAll->Write();
861 d->Write();
862 f.Write();
863 f.Close();
864
865 if (fDiagStep1) fDiagStep1->Reset();
866 if (fDiagStep2) fDiagStep2->Reset();
867 if (fDiagStep3) fDiagStep3->Reset();
868 if (fDiagStep4) fDiagStep4->Reset();
869 if (fDiagAll) fDiagAll->Reset();
121a60bd 870}
871
ddaa8027 872//____________________________________________________________________
873void
874AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
875 AliESDEvent* esd) const
876{
877 TTree* dummy = 0;
878 FillESD(dummy, clusterTree, esd);
879}
cf12b007 880//_____________________________________________________________________
881Bool_t AliFMDReconstructor::GetFMDAsideBit(AliESDFMD* fmd) {
882
883 AliFMDParameters* pars = AliFMDParameters::Instance();
884 Float_t totalMult = 0;
885 for(UShort_t det=1;det<=2;det++) {
886 Int_t nRings = (det == 1 ? 1 : 2);
887 for (UShort_t ir = 0; ir < nRings; ir++) {
888 Char_t ring = (ir == 0 ? 'I' : 'O');
889 UShort_t nsec = (ir == 0 ? 20 : 40);
890 UShort_t nstr = (ir == 0 ? 512 : 256);
891 for(UShort_t sec =0; sec < nsec; sec++) {
892 for(UShort_t strip = 0; strip < nstr; strip++) {
893 Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
894 if(mult == AliESDFMD::kInvalidMult) continue;
895
896 if(mult > pars->GetOfflineTriggerLowCut())
897 totalMult = totalMult + mult;
898 else
899 {
900 if( totalMult > pars->GetOfflineTriggerHitCut()) {
901 return kTRUE;
902 }
903 else totalMult = 0 ;
904 }
905 }
906 }
907 }
908 }
909 return kFALSE;
910
911}
912//_____________________________________________________________________
913Bool_t AliFMDReconstructor::GetFMDCsideBit(AliESDFMD* fmd) {
914
915 AliFMDParameters* pars = AliFMDParameters::Instance();
916 Float_t totalMult = 0;
917 UShort_t det = 3;
918 Int_t nRings = 2;
919 for (UShort_t ir = 0; ir < nRings; ir++) {
920 Char_t ring = (ir == 0 ? 'I' : 'O');
921 UShort_t nsec = (ir == 0 ? 20 : 40);
922 UShort_t nstr = (ir == 0 ? 512 : 256);
923 for(UShort_t sec =0; sec < nsec; sec++) {
924 for(UShort_t strip = 0; strip < nstr; strip++) {
925 Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
926 if(mult == AliESDFMD::kInvalidMult) continue;
927
928 if(mult > pars->GetOfflineTriggerLowCut())
929 totalMult = totalMult + mult;
930 else
931 {
932 if( totalMult > pars->GetOfflineTriggerHitCut()) {
933 return kTRUE;
934 }
935 else totalMult = 0 ;
936 }
937 }
938 }
939 }
940
941 return kFALSE;
942}
42403906 943//____________________________________________________________________
944//
945// EOF
946//