Coding conventions (Christian)
[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 **************************************************************************/
0d0e6995 15/* $Id$ */
c2fc1258 16/** @file AliFMDReconstructor.cxx
17 @author Christian Holm Christensen <cholm@nbi.dk>
18 @date Mon Mar 27 12:47:09 2006
19 @brief FMD reconstruction
20*/
4347b38f 21//____________________________________________________________________
42403906 22//
02a27b50 23// This is a class that constructs AliFMDRecPoint objects from of Digits
4347b38f 24//
25// This class reads either digits from a TClonesArray or raw data from
26// a DDL file (or similar), and stores the read ADC counts in an
27// internal cache (fAdcs).
28//
37c55dc0 29//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
4347b38f 30// Latest changes by Christian Holm Christensen <cholm@nbi.dk>
31//
32//
33//____________________________________________________________________
37c55dc0 34
56b1929b 35#include <AliLog.h> // ALILOG_H
36#include <AliRun.h> // ALIRUN_H
37#include <AliRunLoader.h> // ALIRUNLOADER_H
56b1929b 38#include <AliHeader.h> // ALIHEADER_H
56b1929b 39#include <AliGenEventHeader.h> // ALIGENEVENTHEADER_H
1a1fdef7 40#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
41#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
56b1929b 42#include "AliFMDDigit.h" // ALIFMDDIGIT_H
43#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
56b1929b 44#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
bf000c32 45#include "AliFMDRecPoint.h" // ALIFMDMULTNAIIVE_H
a3537838 46#include "AliESD.h" // ALIESD_H
8f6ee336 47#include <AliESDFMD.h> // ALIESDFMD_H
02a27b50 48class AliRawReader;
4347b38f 49
50//____________________________________________________________________
925e6570 51ClassImp(AliFMDReconstructor)
1a1fdef7 52#if 0
53 ; // This is here to keep Emacs for indenting the next line
54#endif
37c55dc0 55
4347b38f 56//____________________________________________________________________
57AliFMDReconstructor::AliFMDReconstructor()
58 : AliReconstructor(),
e9fd1e20 59 fMult(0x0),
60 fNMult(0),
61 fTreeR(0x0),
62 fCurrentVertex(0),
63 fESDObj(0x0),
64 fESD(0x0)
4347b38f 65{
8f6ee336 66 // Make a new FMD reconstructor object - default CTOR.
4347b38f 67}
68
42403906 69
70//____________________________________________________________________
71AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
8f6ee336 72 : AliReconstructor(),
73 fMult(other.fMult),
e9fd1e20 74 fNMult(other.fNMult),
75 fTreeR(other.fTreeR),
76 fCurrentVertex(other.fCurrentVertex),
77 fESDObj(other.fESDObj),
78 fESD(other.fESD)
42403906 79{
56b1929b 80 // Copy constructor
42403906 81}
82
83
84//____________________________________________________________________
85AliFMDReconstructor&
86AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
87{
56b1929b 88 // Assignment operator
8f6ee336 89 fMult = other.fMult;
e9fd1e20 90 fNMult = other.fNMult;
91 fTreeR = other.fTreeR;
92 fCurrentVertex = other.fCurrentVertex;
8f6ee336 93 fESDObj = other.fESDObj;
e9fd1e20 94 fESD = other.fESD;
42403906 95 return *this;
96}
56b1929b 97
98//____________________________________________________________________
99AliFMDReconstructor::~AliFMDReconstructor()
100{
101 // Destructor
8f6ee336 102 if (fMult) fMult->Delete();
103 if (fMult) delete fMult;
104 if (fESDObj) delete fESDObj;
56b1929b 105}
4347b38f 106
107//____________________________________________________________________
108void
1a1fdef7 109AliFMDReconstructor::Init(AliRunLoader* runLoader)
110{
111 // Initialize the reconstructor
112 AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
bf000c32 113 // Initialize the geometry
9b48326f 114 AliFMDGeometry* geom = AliFMDGeometry::Instance();
115 geom->Init();
116 geom->InitTransformations();
117
118 // Initialize the parameters
119 AliFMDParameters* param = AliFMDParameters::Instance();
120 param->Init();
bf000c32 121
8f6ee336 122 // Current vertex position
1a1fdef7 123 fCurrentVertex = 0;
8f6ee336 124 // Create array of reconstructed strip multiplicities
bf000c32 125 fMult = new TClonesArray("AliFMDRecPoint", 51200);
8f6ee336 126 // Create ESD output object
127 fESDObj = new AliESDFMD;
128
129 // Check that we have a run loader
1a1fdef7 130 if (!runLoader) {
131 Warning("Init", "No run loader");
4347b38f 132 return;
133 }
8f6ee336 134
135 // Check if we can get the header tree
1a1fdef7 136 AliHeader* header = runLoader->GetHeader();
137 if (!header) {
138 Warning("Init", "No header");
4347b38f 139 return;
140 }
8f6ee336 141
142 // Check if we can get a simulation header
1a1fdef7 143 AliGenEventHeader* eventHeader = header->GenEventHeader();
a3537838 144 if (eventHeader) {
145 TArrayF vtx;
146 eventHeader->PrimaryVertex(vtx);
147 fCurrentVertex = vtx[2];
148 AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f",
149 header->GetRun(), header->GetEvent(), fCurrentVertex));
150 Warning("Init", "no generator event header");
151 }
152 else {
153 Warning("Init", "No generator event header - "
154 "perhaps we get the vertex from ESD?");
4347b38f 155 }
4347b38f 156}
dc8af42e 157
4347b38f 158//____________________________________________________________________
159void
1a1fdef7 160AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
161 TTree* digitsTree) const
162{
163 // Convert Raw digits to AliFMDDigit's in a tree
8f6ee336 164 AliDebug(1, "Reading raw data into digits tree");
1a1fdef7 165 AliFMDRawReader rawRead(reader, digitsTree);
166 // rawRead.SetSampleRate(fFMD->GetSampleRate());
167 rawRead.Exec();
4347b38f 168}
169
4347b38f 170//____________________________________________________________________
171void
1a1fdef7 172AliFMDReconstructor::Reconstruct(TTree* digitsTree,
173 TTree* clusterTree) const
4347b38f 174{
1a1fdef7 175 // Reconstruct event from digits in tree
4347b38f 176 // Get the FMD branch holding the digits.
69b696b9 177 // FIXME: The vertex may not be known yet, so we may have to move
178 // some of this to FillESD.
1a1fdef7 179 AliDebug(1, "Reconstructing from digits in a tree");
1e8f773e 180#if 1
a3537838 181 if (fESD) {
8f6ee336 182 const AliESDVertex* vertex = fESD->GetVertex();
183 if (vertex) {
184 AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
185 fCurrentVertex = vertex->GetZv();
186 }
a3537838 187 }
54e415a8 188#endif
1a1fdef7 189 TBranch *digitBranch = digitsTree->GetBranch("FMD");
4347b38f 190 if (!digitBranch) {
1a1fdef7 191 Error("Exec", "No digit branch for the FMD found");
192 return;
4347b38f 193 }
0abe7182 194 TClonesArray* digits = new TClonesArray("AliFMDDigit");
1a1fdef7 195 digitBranch->SetAddress(&digits);
4347b38f 196
8f6ee336 197 if (fMult) fMult->Clear();
198 if (fESDObj) fESDObj->Clear();
199
200 fNMult = 0;
201 fTreeR = clusterTree;
202 fTreeR->Branch("FMD", &fMult);
203
204 AliDebug(5, "Getting entry 0 from digit branch");
1a1fdef7 205 digitBranch->GetEntry(0);
206
8f6ee336 207 AliDebug(5, "Processing digits");
e802be3e 208 ProcessDigits(digits);
209
8f6ee336 210 Int_t written = clusterTree->Fill();
211 AliDebug(10, Form("Filled %d bytes into cluster tree", written));
0abe7182 212 digits->Delete();
213 delete digits;
4347b38f 214}
1a1fdef7 215
8f6ee336 216
4347b38f 217//____________________________________________________________________
e802be3e 218void
219AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f 220{
69b696b9 221 // For each digit, find the pseudo rapdity, azimuthal angle, and
222 // number of corrected ADC counts, and pass it on to the algorithms
223 // used.
e802be3e 224 Int_t nDigits = digits->GetEntries();
1a1fdef7 225 AliDebug(1, Form("Got %d digits", nDigits));
e802be3e 226 for (Int_t i = 0; i < nDigits; i++) {
227 AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
8f6ee336 228 AliFMDParameters* param = AliFMDParameters::Instance();
229 // Check that the strip is not marked as dead
230 if (param->IsDead(digit->Detector(), digit->Ring(),
231 digit->Sector(), digit->Strip())) continue;
232
233 // digit->Print();
234 // Get eta and phi
235 Float_t eta, phi;
236 PhysicalCoordinates(digit, eta, phi);
4347b38f 237
8f6ee336 238 // Substract pedestal.
e802be3e 239 UShort_t counts = SubtractPedestal(digit);
4347b38f 240
8f6ee336 241 // Gain match digits.
242 Double_t edep = Adc2Energy(digit, eta, counts);
243
244 // Make rough multiplicity
245 Double_t mult = Energy2Multiplicity(digit, edep);
246
247 AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
248 "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
249 digit->Detector(), digit->Ring(), digit->Sector(),
250 digit->Strip(), digit->Counts(), counts, edep, mult));
251
252 // Create a `RecPoint' on the output branch.
bf000c32 253 AliFMDRecPoint* m =
254 new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(),
255 digit->Ring(),
256 digit->Sector(),
257 digit->Strip(),
258 eta, phi,
259 edep, mult);
8f6ee336 260 (void)m; // Suppress warnings about unused variables.
261 fNMult++;
262
263 fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(),
264 digit->Sector(), digit->Strip(), mult);
265 fESDObj->SetEta(digit->Detector(), digit->Ring(),
266 digit->Sector(), digit->Strip(), eta);
4347b38f 267 }
4347b38f 268}
8f6ee336 269
1a1fdef7 270//____________________________________________________________________
271UShort_t
272AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
273{
274 // Member function to subtract the pedestal from a digit
275 // This implementation does nothing, but a derived class could over
276 // load this to subtract a pedestal that was given in a database or
277 // something like that.
278
8f6ee336 279 Int_t counts = 0;
280 AliFMDParameters* param = AliFMDParameters::Instance();
281 Float_t pedM = param->GetPedestal(digit->Detector(),
282 digit->Ring(),
283 digit->Sector(),
284 digit->Strip());
285 AliDebug(10, Form("Subtracting pedestal %f from signal %d",
286 pedM, digit->Counts()));
1a1fdef7 287 if (digit->Count3() > 0) counts = digit->Count3();
288 else if (digit->Count2() > 0) counts = digit->Count2();
289 else counts = digit->Count1();
8f6ee336 290 counts = TMath::Max(Int_t(counts - pedM), 0);
291 if (counts > 0) AliDebug(10, "Got a hit strip");
292
1a1fdef7 293 return UShort_t(counts);
294}
295
4347b38f 296//____________________________________________________________________
8f6ee336 297Float_t
298AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit,
299 Float_t /* eta */,
300 UShort_t count) const
301{
302 // Converts number of ADC counts to energy deposited.
303 // Note, that this member function can be overloaded by derived
304 // classes to do strip-specific look-ups in databases or the like,
305 // to find the proper gain for a strip.
306 //
307 // In this simple version, we calculate the energy deposited as
308 //
309 // EnergyDeposited = cos(theta) * gain * count
310 //
311 // where
312 //
313 // Pre_amp_MIP_Range
314 // gain = ----------------- * Energy_deposited_per_MIP
315 // ADC_channel_size
316 //
317 // is constant and the same for all strips.
318
319 // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
320 // Double_t edep = TMath::Abs(TMath::Cos(theta)) * fGain * count;
321 AliFMDParameters* param = AliFMDParameters::Instance();
322 Float_t gain = param->GetPulseGain(digit->Detector(),
323 digit->Ring(),
324 digit->Sector(),
325 digit->Strip());
326 Double_t edep = count * gain;
327 AliDebug(10, Form("Converting counts %d to energy via factor %f",
328 count, gain));
329 return edep;
330}
331
332//____________________________________________________________________
333Float_t
334AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */,
335 Float_t edep) const
336{
337 // Converts an energy signal to number of particles.
338 // Note, that this member function can be overloaded by derived
339 // classes to do strip-specific look-ups in databases or the like,
340 // to find the proper gain for a strip.
341 //
342 // In this simple version, we calculate the multiplicity as
343 //
344 // multiplicity = Energy_deposited / Energy_deposited_per_MIP
345 //
346 // where
347 //
348 // Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
349 //
350 // is constant and the same for all strips
351 AliFMDParameters* param = AliFMDParameters::Instance();
352 Double_t edepMIP = param->GetEdepMip();
353 Float_t mult = edep / edepMIP;
354 if (edep > 0)
355 AliDebug(10, Form("Translating energy %f to multiplicity via "
356 "divider %f->%f", edep, edepMIP, mult));
357 return mult;
358}
359
360//____________________________________________________________________
361void
362AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit,
363 Float_t& eta,
364 Float_t& phi) const
365{
366 // Get the eta and phi of a digit
367 //
368 // Get geometry.
9b48326f 369 AliFMDGeometry* geom = AliFMDGeometry::Instance();
bf000c32 370 Double_t x, y, z, r, theta;
9b48326f 371 geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(),
bf000c32 372 digit->Strip(), x, y, z);
373 // Correct for vertex offset.
374 z += fCurrentVertex;
375 phi = TMath::ATan2(y, x);
376 r = TMath::Sqrt(y * y + x * x);
377 theta = TMath::ATan2(r, z);
378 eta = -TMath::Log(TMath::Tan(theta / 2));
8f6ee336 379}
380
381
382
383//____________________________________________________________________
4347b38f 384void
1a1fdef7 385AliFMDReconstructor::FillESD(TTree* /* digitsTree */,
386 TTree* /* clusterTree */,
8f6ee336 387 AliESD* esd) const
121a60bd 388{
42403906 389 // nothing to be done
69b696b9 390 // FIXME: The vertex may not be known when Reconstruct is executed,
391 // so we may have to move some of that member function here.
8f6ee336 392 AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
393 // fESDObj->Print();
394
395 if (esd) {
396 AliDebug(1, Form("Writing FMD data to ESD tree"));
397 esd->SetFMDData(fESDObj);
a3537838 398 }
121a60bd 399}
400
54e415a8 401
402//____________________________________________________________________
403void
404AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const
405{
406 // Cannot be used. See member function with same name but with 2
407 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 408 AliDebug(1, Form("Calling FillESD with loader and tree"));
54e415a8 409 AliError("MayNotUse");
410}
411//____________________________________________________________________
412void
413AliFMDReconstructor::Reconstruct(AliRunLoader*) const
414{
415 // Cannot be used. See member function with same name but with 2
416 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 417 AliDebug(1, Form("Calling FillESD with loader"));
54e415a8 418 AliError("MayNotUse");
419}
420//____________________________________________________________________
421void
422AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const
423{
424 // Cannot be used. See member function with same name but with 2
425 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 426 AliDebug(1, Form("Calling FillESD with loader and raw reader"));
54e415a8 427 AliError("MayNotUse");
428}
429//____________________________________________________________________
430void
431AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const
432{
433 // Cannot be used. See member function with same name but with 2
434 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 435 AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
54e415a8 436 AliError("MayNotUse");
437}
438//____________________________________________________________________
439void
440AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
441{
442 // Cannot be used. See member function with same name but with 2
443 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 444 AliDebug(1, Form("Calling FillESD with loader and ESD"));
54e415a8 445 AliError("MayNotUse");
446}
447//____________________________________________________________________
448void
449AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const
450{
451 // Cannot be used. See member function with same name but with 2
452 // TTree arguments. Make sure you do local reconstrucion
8f6ee336 453 AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
54e415a8 454 AliError("MayNotUse");
455}
456
42403906 457//____________________________________________________________________
458//
459// EOF
460//