New geometry: SDD, cables and update on V11 (L. Gaudichet)
[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
e802be3e 46#include "AliFMD.h" // ALIFMD_H
47#include "AliFMDDigit.h" // ALIFMDDIGIT_H
48#include "AliFMDParticles.h" // ALIFMDPARTICLES_H
49#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
50#include "AliAltroBuffer.h" // ALIALTROBUFFER_H
51#include "AliLog.h" // ALILOG_H
52#include "AliRun.h" // ALIRUN_H
53#include "AliRunLoader.h" // ALIRUNLOADER_H
54#include "AliLoader.h" // ALILOADER_H
55#include "AliHeader.h" // ALIHEADER_H
56#include "AliGenEventHeader.h" // ALIGENEVENTHEADER_H
57#include "AliFMDRawStream.h" // ALIFMDRAWSTREAM_H
58#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
59#include "AliRawReader.h" // ALIRAWREADER_H
60#include "AliFMDReconstructionAlgorithm.h" // ALIFMDRECONSTRUCTIONALGORITHM_H
4347b38f 61
62//____________________________________________________________________
63ClassImp(AliFMDReconstructor);
37c55dc0 64
4347b38f 65//____________________________________________________________________
66AliFMDReconstructor::AliFMDReconstructor()
67 : AliReconstructor(),
4347b38f 68 fDeltaEta(0),
69 fDeltaPhi(0),
70 fThreshold(0),
71 fPedestal(0),
e802be3e 72 fPedestalWidth(0),
73 fPedestalFactor(0)
4347b38f 74{
42403906 75 // Make a new FMD reconstructor object - default CTOR.
4347b38f 76 SetDeltaEta();
77 SetDeltaPhi();
78 SetThreshold();
79 SetPedestal();
80
81 fParticles = new TClonesArray("AliFMDParticles", 1000);
82 fFMDLoader = 0;
83 fRunLoader = 0;
84 fFMD = 0;
85}
86
42403906 87
88//____________________________________________________________________
89AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
90 : AliReconstructor(),
42403906 91 fDeltaEta(0),
92 fDeltaPhi(0),
93 fThreshold(0),
94 fPedestal(0),
e802be3e 95 fPedestalWidth(0),
96 fPedestalFactor(0)
42403906 97{
98 // Make a new FMD reconstructor object - default CTOR.
99 SetDeltaEta(other.fDeltaEta);
100 SetDeltaPhi(other.fDeltaPhi);
101 SetThreshold(other.fThreshold);
e802be3e 102 SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
42403906 103
104 // fParticles = new TClonesArray("AliFMDParticles", 1000);
105 fFMDLoader = other.fFMDLoader;
106 fRunLoader = other.fRunLoader;
107 fFMD = other.fFMD;
108}
109
110
111//____________________________________________________________________
112AliFMDReconstructor&
113AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
114{
115 // Make a new FMD reconstructor object - default CTOR.
116 SetDeltaEta(other.fDeltaEta);
117 SetDeltaPhi(other.fDeltaPhi);
118 SetThreshold(other.fThreshold);
119 SetPedestal(other.fPedestal, other.fPedestalWidth);
120
121 // fParticles = new TClonesArray("AliFMDParticles", 1000);
122 fFMDLoader = other.fFMDLoader;
123 fRunLoader = other.fRunLoader;
124 fFMD = other.fFMD;
125
126 return *this;
127}
128
4347b38f 129//____________________________________________________________________
130void
e802be3e 131AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor)
4347b38f 132{
133 // Set the pedestal, and pedestal width
e802be3e 134 fPedestal = mean;
135 fPedestalWidth = width;
136 fPedestalFactor = factor;
4347b38f 137}
138
139//____________________________________________________________________
140void
141AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
142 AliRawReader* rawReader) const
143{
144 // Collects all digits in the same active volume into number of
145 // particles
146 //
147 // Reconstruct number of particles in given group of pads for given
148 // FMDvolume determined by numberOfVolume,
149 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
150 // numberOgMaxRing
151 //
152 // The reconstruction method is choosen based on the number of empty
153 // strips.
154 fParticles->Clear();
155 if (!runLoader) {
156 Error("Exec","Run Loader loader is NULL - Session not opened");
157 return;
158 }
159 fRunLoader = runLoader;
160 fFMDLoader = runLoader->GetLoader("FMDLoader");
161 if (!fFMDLoader)
162 Fatal("AliFMDReconstructor","Can not find FMD (loader) "
163 "in specified event");
dc8af42e 164
4347b38f 165 // Get the AliRun object
166 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
167
168 // Get the AliFMD object
169 fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
170 if (!fFMD) {
171 AliError("Can not get FMD from gAlice");
172 return;
173 }
174 fFMDLoader->LoadRecPoints("RECREATE");
dc8af42e 175
4347b38f 176 if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
88cb7938 177
4347b38f 178 if (rawReader) {
179 Int_t event = 0;
180 while (rawReader->NextEvent()) {
e802be3e 181 ProcessEvent(event, rawReader);
4347b38f 182 event++;
183 }
184 }
185 else {
186 Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
187 for(Int_t event = 0; event < nEvents; event++)
e802be3e 188 ProcessEvent(event, 0);
4347b38f 189 }
88cb7938 190
dc8af42e 191
4347b38f 192 fFMDLoader->UnloadRecPoints();
193 fFMDLoader = 0;
194 fRunLoader = 0;
195 fFMD = 0;
196}
dc8af42e 197
4347b38f 198//____________________________________________________________________
199void
200AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e 201{
4347b38f 202 // Collects all digits in the same active volume into number of
203 // particles
204 //
205 // Reconstruct number of particles in given group of pads for given
206 // FMDvolume determined by numberOfVolume,
207 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
208 // numberOgMaxRing
209 //
210 // The reconstruction method is choosen based on the number of empty
211 // strips.
212 Reconstruct(runLoader, 0);
213}
214
215
216//____________________________________________________________________
217void
218AliFMDReconstructor::ProcessEvent(Int_t event,
e802be3e 219 AliRawReader* reader) const
4347b38f 220{
42403906 221 // Process one event read from either a clones array or from a a raw
222 // data reader.
4347b38f 223 fRunLoader->GetEvent(event) ;
224 //event z-vertex for correction eta-rad dependence
225 AliHeader *header = fRunLoader->GetHeader();
e802be3e 226 if (!header) Warning("ProcessEvent", "no AliHeader found!");
4347b38f 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);
e802be3e 232 else Warning("ProcessEvent", "No GenEventHeader Found");
233 fCurrentVertex = o.At(2);
cb1df35e 234
4347b38f 235 // If the recontruction tree isn't loaded, load it
236 if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
88cb7938 237
4347b38f 238 //Make branches to hold the reconstructed particles
239 const Int_t kBufferSize = 16000;
240 fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
241
242 // Load or recreate the digits
243 if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
244 if (!reader) {
245 Error("Exec","Error occured while loading digits. Exiting.");
246 return;
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");
e802be3e 263 TClonesArray* digits = fFMD->Digits();
4347b38f 264 if (!digitBranch) {
265 if (!reader) {
266 Error("Exec", "No digit branch for the FMD found");
267 return;
268 }
269 fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
270 }
271 if (!reader) digitBranch->SetAddress(&digits);
272
e802be3e 273 if (reader) {
274 AliFMDRawReader rawRead(fFMD, reader);
275 // rawRead->SetSampleRate(fSampleRate);
276 rawRead.Exec();
277 }
278 else {
279 // Read the ADC values from a clones array.
280 AliDebug(10, "Reading ADCs from Digits array");
281 // read Digits, and reconstruct the particles
282 if (!fFMDLoader->TreeD()->GetEvent(0)) return;
283 }
4347b38f 284
e802be3e 285 TIter next(&fAlgorithms);
286 AliFMDReconstructionAlgorithm* algorithm = 0;
287 while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
288 algorithm->PreEvent();
4347b38f 289
e802be3e 290 ProcessDigits(digits);
291
292 next.Reset();
293 algorithm = 0;
294 while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
295 algorithm->PostEvent();
296
4347b38f 297 if (reader) {
298 digitTree->Fill();
299 fFMDLoader->WriteDigits("OVERWRITE");
300 }
301 fFMDLoader->UnloadDigits();
302 fFMDLoader->TreeR()->Reset();
303 fFMDLoader->TreeR()->Fill();
304 fFMDLoader->WriteRecPoints("OVERWRITE");
305}
306
307//____________________________________________________________________
e802be3e 308UShort_t
309AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
4347b38f 310{
e802be3e 311 // Member function to subtract the pedestal from a digit
312 // This implementation does nothing, but a derived class could over
313 // load this to subtract a pedestal that was given in a database or
314 // something like that.
4347b38f 315
e802be3e 316 Int_t counts = 0;
317 Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
318 if (digit->Count3() >= 0) counts = digit->Count3();
319 else if (digit->Count2() >= 0) counts = digit->Count2();
320 else counts = digit->Count2();
321 counts = TMath::Max(Int_t(counts - ped), 0);
322 return UShort_t(counts);
4347b38f 323}
324
325//____________________________________________________________________
e802be3e 326void
327AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f 328{
e802be3e 329 Int_t nDigits = digits->GetEntries();
330 for (Int_t i = 0; i < nDigits; i++) {
331 AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
332 AliFMDSubDetector* subDetector = 0;
333 switch (digit->Detector()) {
334 case 1: subDetector = fFMD->GetFMD1(); break;
335 case 2: subDetector = fFMD->GetFMD2(); break;
336 case 3: subDetector = fFMD->GetFMD3(); break;
337 }
338 if (!subDetector) {
339 Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
340 continue;
341 }
4347b38f 342
e802be3e 343 AliFMDRing* ring = 0;
344 Float_t ringZ = 0;
345 switch(digit->Ring()) {
346 case 'i':
347 case 'I':
348 ring = subDetector->GetInner();
349 ringZ = subDetector->GetInnerZ();
350 break;
351 case 'o':
352 case 'O':
353 ring = subDetector->GetOuter();
354 ringZ = subDetector->GetOuterZ();
355 break;
356 }
357 if (!ring) {
358 Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
359 digit->Ring());
360 break;
4347b38f 361 }
362
e802be3e 363 Float_t realZ = fCurrentVertex + ringZ;
364 Float_t stripR = ((ring->GetHighR() - ring->GetLowR()) / ring->GetNStrips()
365 * (digit->Strip() + .5) + ring->GetLowR());
366 Float_t theta = TMath::ATan2(stripR, realZ);
367 Float_t phi = (2 * TMath::Pi() / ring->GetNSectors()
368 * (digit->Sector() + .5));
369 Float_t eta = -TMath::Log(TMath::Tan(theta / 2));
370 UShort_t counts = SubtractPedestal(digit);
4347b38f 371
e802be3e 372 TIter next(&fAlgorithms);
373 AliFMDReconstructionAlgorithm* algorithm = 0;
374 while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
375 algorithm->ProcessDigit(digit, eta, phi, counts);
4347b38f 376 }
4347b38f 377}
1a42d4f2 378
4347b38f 379
4347b38f 380//____________________________________________________________________
381void
e802be3e 382AliFMDReconstructor::ReconstructFromCache() const
4347b38f 383{
384 // Based on the information in the cache, do the reconstruction.
385 Int_t nRecon = 0;
386 // Loop over the detectors
387 for (Int_t i = 1; i <= 3; i++) {
388 AliFMDSubDetector* sub = 0;
389 switch (i) {
390 case 1: sub = fFMD->GetFMD1(); break;
391 case 2: sub = fFMD->GetFMD2(); break;
392 case 3: sub = fFMD->GetFMD3(); break;
393 }
394 if (!sub) continue;
395
396 // Loop over the rings in the detector
397 for (Int_t j = 0; j < 2; j++) {
398 Float_t rZ = 0;
399 AliFMDRing* r = 0;
400 switch (j) {
401 case 0: r = sub->GetInner(); rZ = sub->GetInnerZ(); break;
402 case 1: r = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
403 }
404 if (!r) continue;
405
406 // Calculate low/high theta and eta
407 // FIXME: Is this right?
e802be3e 408 Float_t realZ = fCurrentVertex + rZ;
4347b38f 409 Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
410 Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
411 Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
412 Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
413 if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
414 Float_t tmp = etaIn;
415 etaIn = etaOut;
416 etaOut = tmp;
417 }
418
419 //-------------------------------------------------------------
420 //
421 // Here starts poisson method
422 //
423 // Calculate eta step per strip, number of eta steps, number of
424 // phi steps, and check the sign of the eta increment
425 Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
426 Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
427 Int_t nPhi = Int_t(360. / fDeltaPhi);
428 Float_t sign = TMath::Sign(Float_t(1.), etaIn);
121a60bd 429
4347b38f 430 AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
431 sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
432
433 // Loop over relevant phi values
434 for (Int_t p = 0; p < nPhi; p++) {
435 Float_t minPhi = p * fDeltaPhi;
436 Float_t maxPhi = minPhi + fDeltaPhi;
437 UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
438 UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
439
440 AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
441 minPhi, maxPhi, minSector, maxSector));
442 // Loop over relevant eta values
443 for (Int_t e = nEta; e >= 0; --e) {
444 Float_t maxEta = etaIn - sign * e * fDeltaEta;
445 Float_t minEta = maxEta - sign * fDeltaEta;
446 if (sign > 0) minEta = TMath::Max(minEta, etaOut);
447 else minEta = TMath::Min(minEta, etaOut);
448 Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
449 Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
450 Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
451 Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
452 UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
453 UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
454
455 AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
456 " [radii %f, %f, thetas %f, %f, sign %d]",
457 minEta, maxEta, minStrip, maxStrip,
458 minR, maxR, theta1, theta2, sign));
459
460 // Count number of empty strips
461 Int_t emptyStrips = 0;
462 for (Int_t sector = minSector; sector < maxSector; sector++)
e802be3e 463 for (Int_t strip = minStrip; strip < maxStrip; strip++) emptyStrips++;
464 // if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip)
465 // < fThreshold) emptyStrips++;
4347b38f 466
467 // The total number of strips
468 Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
469
470 // Log ratio of empty to total number of strips
471 AliDebug(10, Form("Lambda= %d / %d = %f",
472 emptyStrips, nTotal,
473 Float_t(emptyStrips) / nTotal));
474
475 Double_t lambda = (emptyStrips > 0 ?
476 - TMath::Log(Double_t(emptyStrips) / nTotal) :
477 1);
478
479 // The reconstructed number of particles is then given by
480 Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
481
482 // Add a AliFMDParticles to the reconstruction tree.
483 new((*fParticles)[nRecon])
484 AliFMDParticles(sub->GetId(), r->GetId(),
485 minSector, maxSector, minStrip, maxStrip,
486 minEta, maxEta, minPhi, maxPhi,
487 reconstructed, AliFMDParticles::kPoission);
488 nRecon++;
489 } // phi
490 } // eta
491 } // ring
492 } // detector
493}
494
495
496//____________________________________________________________________
497void
498AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
499 AliESD* /*esd*/) const
121a60bd 500{
42403906 501 // nothing to be done
121a60bd 502
503}
504
42403906 505//____________________________________________________________________
506//
507// EOF
508//