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