Implement the new logging system (AliLog)
[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 **************************************************************************/
15
4347b38f 16//____________________________________________________________________
17// This is a class that constructs ReconstParticles (reconstructed
18// particles) out of Digits
19//
20//
21// This class reads either digits from a TClonesArray or raw data from
22// a DDL file (or similar), and stores the read ADC counts in an
23// internal cache (fAdcs).
24//
25// From the cached values it then calculates the number of particles
26// that hit a region of the FMDs, as specified by the user.
27//
28// The reconstruction can be done in two ways: Either via counting the
29// number of empty strips (Poisson method), or by converting the ADC
30// signal to an energy deposition, and then dividing by the typical
31// energy loss of a particle.
32//
33// Currently, this class only reads the digits from a TClonesArray,
34// and the Poission method for reconstruction.
35//
37c55dc0 36//
37//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f 38// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
39//
40//
41//____________________________________________________________________
37c55dc0 42
8f1cfb0c 43#include "AliFMD.h"
44#include "AliFMDDigit.h"
45#include "AliFMDParticles.h"
46#include "AliFMDReconstructor.h"
47#include "AliAltroBuffer.h"
48#include "AliLog.h"
49#include "AliRun.h"
50#include "AliRunLoader.h"
51#include "AliLoader.h"
52#include "AliHeader.h"
53#include "AliGenEventHeader.h"
54#include "AliFMDRawStream.h"
55#include "AliRawReader.h"
4347b38f 56
57//____________________________________________________________________
58ClassImp(AliFMDReconstructor);
37c55dc0 59
4347b38f 60//____________________________________________________________________
61AliFMDReconstructor::AliFMDReconstructor()
62 : AliReconstructor(),
63 fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
64 fDeltaEta(0),
65 fDeltaPhi(0),
66 fThreshold(0),
67 fPedestal(0),
68 fPedestalWidth(0)
69{
70 SetDeltaEta();
71 SetDeltaPhi();
72 SetThreshold();
73 SetPedestal();
74
75 fParticles = new TClonesArray("AliFMDParticles", 1000);
76 fFMDLoader = 0;
77 fRunLoader = 0;
78 fFMD = 0;
79}
80
81//____________________________________________________________________
82void
83AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width)
84{
85 // Set the pedestal, and pedestal width
86 fPedestal = mean;
87 fPedestalWidth = width;
88}
89
90//____________________________________________________________________
91void
92AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
93 AliRawReader* rawReader) const
94{
95 // Collects all digits in the same active volume into number of
96 // particles
97 //
98 // Reconstruct number of particles in given group of pads for given
99 // FMDvolume determined by numberOfVolume,
100 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
101 // numberOgMaxRing
102 //
103 // The reconstruction method is choosen based on the number of empty
104 // strips.
105 fParticles->Clear();
106 if (!runLoader) {
107 Error("Exec","Run Loader loader is NULL - Session not opened");
108 return;
109 }
110 fRunLoader = runLoader;
111 fFMDLoader = runLoader->GetLoader("FMDLoader");
112 if (!fFMDLoader)
113 Fatal("AliFMDReconstructor","Can not find FMD (loader) "
114 "in specified event");
dc8af42e 115
4347b38f 116 // Get the AliRun object
117 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
118
119 // Get the AliFMD object
120 fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
121 if (!fFMD) {
122 AliError("Can not get FMD from gAlice");
123 return;
124 }
125 fFMDLoader->LoadRecPoints("RECREATE");
dc8af42e 126
4347b38f 127 if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
88cb7938 128
4347b38f 129 TClonesArray* digits = fFMD->Digits();
130 if (rawReader) {
131 Int_t event = 0;
132 while (rawReader->NextEvent()) {
133 ProcessEvent(event, rawReader, digits);
134 event++;
135 }
136 }
137 else {
138 Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
139 for(Int_t event = 0; event < nEvents; event++)
140 ProcessEvent(event, 0, digits);
141 }
88cb7938 142
dc8af42e 143
4347b38f 144 fFMDLoader->UnloadRecPoints();
145 fFMDLoader = 0;
146 fRunLoader = 0;
147 fFMD = 0;
148}
dc8af42e 149
4347b38f 150//____________________________________________________________________
151void
152AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e 153{
4347b38f 154 // Collects all digits in the same active volume into number of
155 // particles
156 //
157 // Reconstruct number of particles in given group of pads for given
158 // FMDvolume determined by numberOfVolume,
159 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
160 // numberOgMaxRing
161 //
162 // The reconstruction method is choosen based on the number of empty
163 // strips.
164 Reconstruct(runLoader, 0);
165}
166
167
168//____________________________________________________________________
169void
170AliFMDReconstructor::ProcessEvent(Int_t event,
171 AliRawReader* reader,
172 TClonesArray* digits) const
173{
174 fRunLoader->GetEvent(event) ;
175 //event z-vertex for correction eta-rad dependence
176 AliHeader *header = fRunLoader->GetHeader();
177 if (!header)
178 Warning("ProcessEvent", "no AliHeader found!");
179 AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
3d44ce66 180
4347b38f 181 // Get the Z--coordinate from the event header
182 TArrayF o(3);
183 if (genHeader) genHeader->PrimaryVertex(o);
184 Float_t zVertex = o.At(2);
cb1df35e 185
4347b38f 186 // If the recontruction tree isn't loaded, load it
187 if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
88cb7938 188
4347b38f 189 //Make branches to hold the reconstructed particles
190 const Int_t kBufferSize = 16000;
191 fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
192
193 // Load or recreate the digits
194 if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
195 if (!reader) {
196 Error("Exec","Error occured while loading digits. Exiting.");
197 return;
198 }
46501dfb 199
4347b38f 200 }
201 // Get the digits tree
202 TTree* digitTree = fFMDLoader->TreeD();
203 if (!digitTree) {
204 if (!reader) {
205 Error("Exec","Can not get Tree with Digits. "
206 "Nothing to reconstruct - Exiting");
207 return;
208 }
209 fFMDLoader->MakeTree("D");
210 digitTree = fFMDLoader->TreeD();
211
212 }
213 // Get the FMD branch holding the digits.
214 TBranch *digitBranch = digitTree->GetBranch("FMD");
215 if (!digitBranch) {
216 if (!reader) {
217 Error("Exec", "No digit branch for the FMD found");
218 return;
219 }
220 fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
221 }
222 if (!reader) digitBranch->SetAddress(&digits);
223
224 fEmptyStrips = 0;
225 fTotalStrips = 0;
226 Bool_t ok = kFALSE;
227 if (reader) ok = ReadAdcs(reader);
228 else if (digits) ok = ReadAdcs(digits);
229 if (!ok) return;
230
231 ReconstructFromCache(zVertex);
232
233 if (reader) {
234 digitTree->Fill();
235 fFMDLoader->WriteDigits("OVERWRITE");
236 }
237 fFMDLoader->UnloadDigits();
238 fFMDLoader->TreeR()->Reset();
239 fFMDLoader->TreeR()->Fill();
240 fFMDLoader->WriteRecPoints("OVERWRITE");
241}
242
243//____________________________________________________________________
244Bool_t
245AliFMDReconstructor::ReadAdcs(TClonesArray* digits) const
246{
247 AliDebug(10, "Reading ADCs from Digits array");
248 // read Digits, and reconstruct the particles
249 if (!fFMDLoader->TreeD()->GetEvent(0)) return kFALSE;
250
251 // Reads the digits from the array, and fills up the cache (fAdcs)
252 fAdcs.Clear();
253 Int_t nDigits = digits->GetEntries();
254 for (Int_t digit = 0; digit < nDigits; digit++) {
255 AliFMDDigit* fmdDigit =
256 static_cast<AliFMDDigit*>(digits->UncheckedAt(digit));
257
258 ProcessDigit(fmdDigit);
259 } //digit loop
260 return kTRUE;
261}
262
263//____________________________________________________________________
264Bool_t
265AliFMDReconstructor::ReadAdcs(AliRawReader* reader) const
266{
267 AliDebug(10, "Reading ADCs from RawReader");
268 // Reads the digits from a RAW data
269 fAdcs.Clear();
270 // reader->Reset();
271
272 if (!reader->ReadHeader()) {
273 Error("ReadAdcs", "Couldn't read header");
274 return kFALSE;
275 }
276
277 // Use AliAltroRawStream to read the ALTRO format. No need to
278 // reinvent the wheel :-)
279 AliFMDRawStream input(reader);
280 // Select FMD DDL's
281 reader->Select(AliFMD::kBaseDDL >> 8);
282
283 Int_t oldDDL = -1;
284 Int_t count = 0;
285 UShort_t detector = 1; // Must be one here
286 UShort_t oldDetector = 0;
287 // Loop over data in file
288 Bool_t next = kTRUE;
289
290 // local Cache
291 TArrayI counts(10);
292 counts.Reset(-1);
293 Int_t offset = 0;
294
295 while (next) {
296 next = input.Next();
297
298
299 count++;
300 Int_t ddl = reader->GetDDLID();
301 if (ddl != oldDDL
302 || input.IsNewStrip()
303 || !next) {
304 // Make a new digit, if we have some data!
305 if (counts[0] >= 0) {
306 // Got a new strip.
307 AliDebug(10, Form("Add a new strip: FMD%d%c[%2d,%3d] "
308 "(current: FMD%d%c[%2d,%3d])",
309 oldDetector, input.PrevRing(),
310 input.PrevSector() , input.PrevStrip(),
311 detector , input.Ring(), input.Sector(),
312 input.Strip()));
313 fFMD->AddDigit(oldDetector,
314 input.PrevRing(),
315 input.PrevSector(),
316 input.PrevStrip(),
317 counts[0], counts[1], counts[2]);
318 AliFMDDigit* digit =
319 static_cast<AliFMDDigit*>(fFMD->Digits()->
320 UncheckedAt(fFMD->GetNdigits()-1));
321 ProcessDigit(digit);
322 }
323
324 if (!next) {
325 AliDebug(10, Form("Read %d channels for FMD%d",
326 count + 1, detector));
327 break;
328 }
329
330
331 // If we got a new DDL, it means we have a new detector.
332 if (ddl != oldDDL) {
333 if (detector != 0)
334 AliDebug(10, Form("Read %d channels for FMD%d",
335 count + 1, detector));
336 // Reset counts, and update the DDL cache
337 count = 0;
338 oldDDL = ddl;
339 // Check that we're processing a FMD detector
340 Int_t detId = reader->GetDetectorID();
341 if (detId != (AliFMD::kBaseDDL >> 8)) {
342 Error("ReadAdcs", "Detector ID %d != %d",
343 detId, (AliFMD::kBaseDDL >> 8));
344 break;
46501dfb 345 }
4347b38f 346 // Figure out what detector we're deling with
347 oldDetector = detector;
348 switch (ddl) {
349 case 0: detector = 1; break;
350 case 1: detector = 2; break;
351 case 2: detector = 3; break;
352 default:
353 Error("ReadAdcs", "Unknown DDL 0x%x for FMD", ddl);
354 return kFALSE;
46501dfb 355 }
4347b38f 356 AliDebug(10, Form("Reading ADCs for 0x%x - That is FMD%d",
357 reader->GetEquipmentId(), detector));
358 }
359 counts.Reset(-1);
360 offset = 0;
361 }
362
363 counts[offset] = input.Count();
364 offset++;
365
366 AliDebug(10, Form("ADC of FMD%d%c[%2d,%3d] += %d",
367 detector, input.Ring(), input.Sector(),
368 input.Strip(), input.Count()));
369 oldDetector = detector;
370 }
371 return kTRUE;
372}
88cb7938 373
4347b38f 374//____________________________________________________________________
375void
376AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
377{
378 // Process a digit. Derived classes can overload this member
379 // function to do stuff to the digit. However, it should write the
380 // ADC count to the internal cache
381 //
382 // fAdcs(detector - 1, ring, sector, strip) = counts;
383 //
384 // In this implementation, we count the number of strips below
385 // threshold. This we do to later choose what kind of
386 // reconstruction algorithm we'd like to use.
387 //
388 UShort_t detector = digit->Detector();
389 Char_t ring = digit->Ring();
390 UShort_t sector = digit->Sector();
391 UShort_t strip = digit->Strip();
392
393 UShort_t counts = SubtractPedestal(digit);
394
395 fAdcs(detector - 1, ring, sector, strip) = counts;
396 if (counts < fThreshold) fEmptyStrips++;
397 fTotalStrips++;
398}
1a42d4f2 399
4347b38f 400//____________________________________________________________________
401UShort_t
402AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
403{
404 // Member function to subtract the pedestal from a digit
405 // This implementation does nothing, but a derived class could over
406 // load this to subtract a pedestal that was given in a database or
407 // something like that.
3d44ce66 408
4347b38f 409 Int_t counts =
410 TMath::Max(Int_t(digit->Count1() - fPedestal - 3 * fPedestalWidth), 0);
411 if (digit->Count2() >= 0)
412 counts +=
413 TMath::Max(Int_t(digit->Count2() - fPedestal - 3 * fPedestalWidth), 0);
414 if (digit->Count3() >= 0)
415 counts +=
416 TMath::Max(Int_t(digit->Count3() - fPedestal - 3 * fPedestalWidth), 0);
417
418 if (counts < 0) counts = 0;
419 return UShort_t(counts);
88cb7938 420}
dc8af42e 421
4347b38f 422//____________________________________________________________________
423void
424AliFMDReconstructor::ReconstructFromCache(Float_t zVertex) const
425{
426 // Based on the information in the cache, do the reconstruction.
427 Int_t nRecon = 0;
428 // Loop over the detectors
429 for (Int_t i = 1; i <= 3; i++) {
430 AliFMDSubDetector* sub = 0;
431 switch (i) {
432 case 1: sub = fFMD->GetFMD1(); break;
433 case 2: sub = fFMD->GetFMD2(); break;
434 case 3: sub = fFMD->GetFMD3(); break;
435 }
436 if (!sub) continue;
437
438 // Loop over the rings in the detector
439 for (Int_t j = 0; j < 2; j++) {
440 Float_t rZ = 0;
441 AliFMDRing* r = 0;
442 switch (j) {
443 case 0: r = sub->GetInner(); rZ = sub->GetInnerZ(); break;
444 case 1: r = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
445 }
446 if (!r) continue;
447
448 // Calculate low/high theta and eta
449 // FIXME: Is this right?
450 Float_t realZ = zVertex + rZ;
451 Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
452 Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
453 Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
454 Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
455 if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
456 Float_t tmp = etaIn;
457 etaIn = etaOut;
458 etaOut = tmp;
459 }
460
461 //-------------------------------------------------------------
462 //
463 // Here starts poisson method
464 //
465 // Calculate eta step per strip, number of eta steps, number of
466 // phi steps, and check the sign of the eta increment
467 Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
468 Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
469 Int_t nPhi = Int_t(360. / fDeltaPhi);
470 Float_t sign = TMath::Sign(Float_t(1.), etaIn);
121a60bd 471
4347b38f 472 AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
473 sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
474
475 // Loop over relevant phi values
476 for (Int_t p = 0; p < nPhi; p++) {
477 Float_t minPhi = p * fDeltaPhi;
478 Float_t maxPhi = minPhi + fDeltaPhi;
479 UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
480 UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
481
482 AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
483 minPhi, maxPhi, minSector, maxSector));
484 // Loop over relevant eta values
485 for (Int_t e = nEta; e >= 0; --e) {
486 Float_t maxEta = etaIn - sign * e * fDeltaEta;
487 Float_t minEta = maxEta - sign * fDeltaEta;
488 if (sign > 0) minEta = TMath::Max(minEta, etaOut);
489 else minEta = TMath::Min(minEta, etaOut);
490 Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
491 Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
492 Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
493 Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
494 UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
495 UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
496
497 AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
498 " [radii %f, %f, thetas %f, %f, sign %d]",
499 minEta, maxEta, minStrip, maxStrip,
500 minR, maxR, theta1, theta2, sign));
501
502 // Count number of empty strips
503 Int_t emptyStrips = 0;
504 for (Int_t sector = minSector; sector < maxSector; sector++)
505 for (Int_t strip = minStrip; strip < maxStrip; strip++)
506 if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip)
507 < fThreshold) emptyStrips++;
508
509 // The total number of strips
510 Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
511
512 // Log ratio of empty to total number of strips
513 AliDebug(10, Form("Lambda= %d / %d = %f",
514 emptyStrips, nTotal,
515 Float_t(emptyStrips) / nTotal));
516
517 Double_t lambda = (emptyStrips > 0 ?
518 - TMath::Log(Double_t(emptyStrips) / nTotal) :
519 1);
520
521 // The reconstructed number of particles is then given by
522 Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
523
524 // Add a AliFMDParticles to the reconstruction tree.
525 new((*fParticles)[nRecon])
526 AliFMDParticles(sub->GetId(), r->GetId(),
527 minSector, maxSector, minStrip, maxStrip,
528 minEta, maxEta, minPhi, maxPhi,
529 reconstructed, AliFMDParticles::kPoission);
530 nRecon++;
531 } // phi
532 } // eta
533 } // ring
534 } // detector
535}
536
537
538//____________________________________________________________________
539void
540AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
541 AliESD* /*esd*/) const
121a60bd 542{
543// nothing to be done
544
545}
546