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