changed default value for fGeVcharge
[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//
7684b53c 20// This is a class that constructs AliFMDMult (reconstructed
21// multiplicity) from of Digits
4347b38f 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//
7684b53c 35// +---------------------+ +---------------------+
36// | AliFMDReconstructor |<>-----| AliFMDMultAlgorithm |
37// +---------------------+ +---------------------+
38// ^
39// |
40// +-----------+---------+
41// | |
42// +-------------------+ +------------------+
43// | AliFMDMultPoisson | | AliFMDMultNaiive |
44// +-------------------+ +------------------+
45//
46// AliFMDReconstructor acts as a manager class. It contains a list of
47// AliFMDMultAlgorithm objects. The call graph looks something like
48//
49//
50// +----------------------+ +----------------------+
51// | :AliFMDReconstructor | | :AliFMDMultAlgorithm |
52// +----------------------+ +----------------------+
53// | |
54// Reconstruct +-+ |
55// ------------>| | PreRun +-+
56// | |------------------------------->| |
57// | | +-+
58// | |-----+ (for each event) |
59// | | | *ProcessEvent |
60// |+-+ | |
61// || |<---+ PreEvent +-+
62// || |------------------------------>| |
63// || | +-+
64// || |-----+ |
65// || | | ProcessDigits |
66// ||+-+ | |
67// ||| |<---+ |
68// ||| | *ProcessDigit(digit) +-+
69// ||| |----------------------------->| |
70// ||| | +-+
71// ||+-+ |
72// || | PostEvent +-+
73// || |------------------------------>| |
74// || | +-+
75// |+-+ |
76// | | PostRun +-+
77// | |------------------------------->| |
78// | | +-+
79// +-+ |
80// | |
81//
4347b38f 82//
37c55dc0 83//
84//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f 85// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
86//
87//
88//____________________________________________________________________
37c55dc0 89
56b1929b 90#include <AliLog.h> // ALILOG_H
91#include <AliRun.h> // ALIRUN_H
92#include <AliRunLoader.h> // ALIRUNLOADER_H
93#include <AliLoader.h> // ALILOADER_H
94#include <AliHeader.h> // ALIHEADER_H
95#include <AliRawReader.h> // ALIRAWREADER_H
96#include <AliGenEventHeader.h> // ALIGENEVENTHEADER_H
97#include "AliFMD.h" // ALIFMD_H
98#include "AliFMDDigit.h" // ALIFMDDIGIT_H
99#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
100#include "AliFMDRawStream.h" // ALIFMDRAWSTREAM_H
101#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
102#include "AliFMDMultAlgorithm.h" // ALIFMDMULTALGORITHM_H
103#include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
104#include "AliFMDMultNaiive.h" // ALIFMDMULTNAIIVE_H
4347b38f 105
106//____________________________________________________________________
925e6570 107ClassImp(AliFMDReconstructor)
37c55dc0 108
4347b38f 109//____________________________________________________________________
110AliFMDReconstructor::AliFMDReconstructor()
111 : AliReconstructor(),
4347b38f 112 fPedestal(0),
e802be3e 113 fPedestalWidth(0),
114 fPedestalFactor(0)
4347b38f 115{
42403906 116 // Make a new FMD reconstructor object - default CTOR.
4347b38f 117 SetPedestal();
118
4347b38f 119 fFMDLoader = 0;
120 fRunLoader = 0;
121 fFMD = 0;
56b1929b 122 fAlgorithms.Add(new AliFMDMultNaiive);
123 fAlgorithms.Add(new AliFMDMultPoisson);
4347b38f 124}
125
42403906 126
127//____________________________________________________________________
128AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
129 : AliReconstructor(),
42403906 130 fPedestal(0),
e802be3e 131 fPedestalWidth(0),
132 fPedestalFactor(0)
42403906 133{
56b1929b 134 // Copy constructor
e802be3e 135 SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
42403906 136
42403906 137 fFMDLoader = other.fFMDLoader;
138 fRunLoader = other.fRunLoader;
139 fFMD = other.fFMD;
56b1929b 140
141 fAlgorithms.Delete();
142 TIter next(&(other.fAlgorithms));
143 AliFMDMultAlgorithm* algorithm = 0;
144 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
145 fAlgorithms.Add(algorithm);
146 fAlgorithms.SetOwner(kFALSE);
42403906 147}
148
149
150//____________________________________________________________________
151AliFMDReconstructor&
152AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
153{
56b1929b 154 // Assignment operator
155 SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
42403906 156
42403906 157 fFMDLoader = other.fFMDLoader;
158 fRunLoader = other.fRunLoader;
159 fFMD = other.fFMD;
160
56b1929b 161 fAlgorithms.Delete();
162 TIter next(&(other.fAlgorithms));
163 AliFMDMultAlgorithm* algorithm = 0;
164 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
165 fAlgorithms.Add(algorithm);
166 fAlgorithms.SetOwner(kFALSE);
167
42403906 168 return *this;
169}
56b1929b 170
171//____________________________________________________________________
172AliFMDReconstructor::~AliFMDReconstructor()
173{
174 // Destructor
175 fAlgorithms.Delete();
176}
42403906 177
4347b38f 178//____________________________________________________________________
179void
e802be3e 180AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor)
4347b38f 181{
182 // Set the pedestal, and pedestal width
e802be3e 183 fPedestal = mean;
184 fPedestalWidth = width;
185 fPedestalFactor = factor;
4347b38f 186}
187
188//____________________________________________________________________
189void
190AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
191 AliRawReader* rawReader) const
192{
193 // Collects all digits in the same active volume into number of
194 // particles
195 //
196 // Reconstruct number of particles in given group of pads for given
197 // FMDvolume determined by numberOfVolume,
198 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
199 // numberOgMaxRing
200 //
201 // The reconstruction method is choosen based on the number of empty
202 // strips.
4347b38f 203 if (!runLoader) {
204 Error("Exec","Run Loader loader is NULL - Session not opened");
205 return;
206 }
207 fRunLoader = runLoader;
208 fFMDLoader = runLoader->GetLoader("FMDLoader");
209 if (!fFMDLoader)
210 Fatal("AliFMDReconstructor","Can not find FMD (loader) "
211 "in specified event");
dc8af42e 212
4347b38f 213 // Get the AliRun object
214 if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
215
216 // Get the AliFMD object
217 fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
218 if (!fFMD) {
219 AliError("Can not get FMD from gAlice");
220 return;
221 }
222 fFMDLoader->LoadRecPoints("RECREATE");
dc8af42e 223
4347b38f 224 if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
88cb7938 225
56b1929b 226 TIter next(&fAlgorithms);
227 AliFMDMultAlgorithm* algorithm = 0;
228 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
229 algorithm->PreRun(fFMD);
230
4347b38f 231 if (rawReader) {
232 Int_t event = 0;
233 while (rawReader->NextEvent()) {
e802be3e 234 ProcessEvent(event, rawReader);
4347b38f 235 event++;
236 }
237 }
238 else {
239 Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
240 for(Int_t event = 0; event < nEvents; event++)
e802be3e 241 ProcessEvent(event, 0);
4347b38f 242 }
88cb7938 243
56b1929b 244 next.Reset();
245 algorithm = 0;
246 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
247 algorithm->PostRun();
dc8af42e 248
4347b38f 249 fFMDLoader->UnloadRecPoints();
250 fFMDLoader = 0;
251 fRunLoader = 0;
252 fFMD = 0;
253}
dc8af42e 254
4347b38f 255//____________________________________________________________________
256void
257AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
dc8af42e 258{
4347b38f 259 // Collects all digits in the same active volume into number of
260 // particles
261 //
262 // Reconstruct number of particles in given group of pads for given
263 // FMDvolume determined by numberOfVolume,
264 // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
265 // numberOgMaxRing
266 //
267 // The reconstruction method is choosen based on the number of empty
268 // strips.
269 Reconstruct(runLoader, 0);
270}
271
272
273//____________________________________________________________________
274void
275AliFMDReconstructor::ProcessEvent(Int_t event,
e802be3e 276 AliRawReader* reader) const
4347b38f 277{
42403906 278 // Process one event read from either a clones array or from a a raw
279 // data reader.
4347b38f 280 fRunLoader->GetEvent(event) ;
281 //event z-vertex for correction eta-rad dependence
282 AliHeader *header = fRunLoader->GetHeader();
e802be3e 283 if (!header) Warning("ProcessEvent", "no AliHeader found!");
4347b38f 284 AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
3d44ce66 285
4347b38f 286 // Get the Z--coordinate from the event header
287 TArrayF o(3);
288 if (genHeader) genHeader->PrimaryVertex(o);
e802be3e 289 else Warning("ProcessEvent", "No GenEventHeader Found");
290 fCurrentVertex = o.At(2);
cb1df35e 291
4347b38f 292 // If the recontruction tree isn't loaded, load it
293 if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
88cb7938 294
4347b38f 295 // Load or recreate the digits
296 if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
297 if (!reader) {
298 Error("Exec","Error occured while loading digits. Exiting.");
299 return;
300 }
4347b38f 301 }
302 // Get the digits tree
303 TTree* digitTree = fFMDLoader->TreeD();
304 if (!digitTree) {
305 if (!reader) {
306 Error("Exec","Can not get Tree with Digits. "
307 "Nothing to reconstruct - Exiting");
308 return;
309 }
310 fFMDLoader->MakeTree("D");
311 digitTree = fFMDLoader->TreeD();
312
313 }
314 // Get the FMD branch holding the digits.
315 TBranch *digitBranch = digitTree->GetBranch("FMD");
e802be3e 316 TClonesArray* digits = fFMD->Digits();
4347b38f 317 if (!digitBranch) {
318 if (!reader) {
319 Error("Exec", "No digit branch for the FMD found");
320 return;
321 }
322 fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
323 }
324 if (!reader) digitBranch->SetAddress(&digits);
325
e802be3e 326 if (reader) {
327 AliFMDRawReader rawRead(fFMD, reader);
7684b53c 328 AliDebug(10, Form("Making raw reader with sample rate: %d",
329 fFMD->GetSampleRate()));
330 rawRead.SetSampleRate(fFMD->GetSampleRate());
e802be3e 331 rawRead.Exec();
332 }
333 else {
334 // Read the ADC values from a clones array.
335 AliDebug(10, "Reading ADCs from Digits array");
336 // read Digits, and reconstruct the particles
337 if (!fFMDLoader->TreeD()->GetEvent(0)) return;
338 }
4347b38f 339
e802be3e 340 TIter next(&fAlgorithms);
56b1929b 341 AliFMDMultAlgorithm* algorithm = 0;
342 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
343 algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);
4347b38f 344
e802be3e 345 ProcessDigits(digits);
346
347 next.Reset();
348 algorithm = 0;
56b1929b 349 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
e802be3e 350 algorithm->PostEvent();
351
4347b38f 352 if (reader) {
353 digitTree->Fill();
354 fFMDLoader->WriteDigits("OVERWRITE");
355 }
356 fFMDLoader->UnloadDigits();
357 fFMDLoader->TreeR()->Reset();
358 fFMDLoader->TreeR()->Fill();
359 fFMDLoader->WriteRecPoints("OVERWRITE");
360}
361
362//____________________________________________________________________
e802be3e 363UShort_t
364AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
4347b38f 365{
e802be3e 366 // Member function to subtract the pedestal from a digit
367 // This implementation does nothing, but a derived class could over
368 // load this to subtract a pedestal that was given in a database or
369 // something like that.
4347b38f 370
e802be3e 371 Int_t counts = 0;
7684b53c 372 Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
373 if (digit->Count3() > 0) counts = digit->Count3();
374 else if (digit->Count2() > 0) counts = digit->Count2();
375 else counts = digit->Count1();
e802be3e 376 counts = TMath::Max(Int_t(counts - ped), 0);
377 return UShort_t(counts);
4347b38f 378}
379
380//____________________________________________________________________
e802be3e 381void
382AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f 383{
e802be3e 384 Int_t nDigits = digits->GetEntries();
385 for (Int_t i = 0; i < nDigits; i++) {
386 AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
387 AliFMDSubDetector* subDetector = 0;
388 switch (digit->Detector()) {
389 case 1: subDetector = fFMD->GetFMD1(); break;
390 case 2: subDetector = fFMD->GetFMD2(); break;
391 case 3: subDetector = fFMD->GetFMD3(); break;
392 }
393 if (!subDetector) {
394 Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
395 continue;
396 }
4347b38f 397
e802be3e 398 AliFMDRing* ring = 0;
399 Float_t ringZ = 0;
400 switch(digit->Ring()) {
401 case 'i':
402 case 'I':
403 ring = subDetector->GetInner();
404 ringZ = subDetector->GetInnerZ();
405 break;
406 case 'o':
407 case 'O':
408 ring = subDetector->GetOuter();
409 ringZ = subDetector->GetOuterZ();
410 break;
411 }
412 if (!ring) {
413 Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
414 digit->Ring());
415 break;
4347b38f 416 }
417
e802be3e 418 Float_t realZ = fCurrentVertex + ringZ;
56b1929b 419 Float_t stripR = ((ring->GetHighR() - ring->GetLowR())
420 / ring->GetNStrips() * (digit->Strip() + .5)
421 + ring->GetLowR());
e802be3e 422 Float_t theta = TMath::ATan2(stripR, realZ);
423 Float_t phi = (2 * TMath::Pi() / ring->GetNSectors()
424 * (digit->Sector() + .5));
425 Float_t eta = -TMath::Log(TMath::Tan(theta / 2));
426 UShort_t counts = SubtractPedestal(digit);
4347b38f 427
e802be3e 428 TIter next(&fAlgorithms);
56b1929b 429 AliFMDMultAlgorithm* algorithm = 0;
430 while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
e802be3e 431 algorithm->ProcessDigit(digit, eta, phi, counts);
4347b38f 432 }
4347b38f 433}
4347b38f 434
4347b38f 435
436//____________________________________________________________________
437void
438AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
439 AliESD* /*esd*/) const
121a60bd 440{
42403906 441 // nothing to be done
121a60bd 442
443}
444
42403906 445//____________________________________________________________________
446//
447// EOF
448//