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