Added ITS Alignment interface to MillePede2 and related supermodule class
[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
6169f936 24// This class reads either digits from a TClonesArray or raw data from
4347b38f 25// a DDL file (or similar), and stores the read ADC counts in an
6169f936 26// internal cache (fAdcs). The rec-points are made via the naiive
27// method.
4347b38f 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
f95a63c4 35// #include <AliLog.h> // ALILOG_H
6169f936 36// #include <AliRun.h> // ALIRUN_H
f95a63c4 37#include "AliFMDDebug.h"
1a1fdef7 38#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
39#include "AliFMDParameters.h" // ALIFMDPARAMETERS_H
5cf05dbb 40#include "AliFMDAltroMapping.h" // ALIFMDALTROMAPPING_H
56b1929b 41#include "AliFMDDigit.h" // ALIFMDDIGIT_H
42#include "AliFMDReconstructor.h" // ALIFMDRECONSTRUCTOR_H
56b1929b 43#include "AliFMDRawReader.h" // ALIFMDRAWREADER_H
bf000c32 44#include "AliFMDRecPoint.h" // ALIFMDMULTNAIIVE_H
af885e0f 45#include "AliESDEvent.h" // ALIESDEVENT_H
aad72f45 46#include "AliESDVertex.h" // ALIESDVERTEX_H
50b9d194 47#include "AliESDTZERO.h" // ALIESDVERTEX_H
8f6ee336 48#include <AliESDFMD.h> // ALIESDFMD_H
aad72f45 49#include <TMath.h>
9684be2f 50#include <TH1.h>
51#include <TH2.h>
52#include <TFile.h>
68aba90a 53#include <climits>
50b9d194 54// Import revertexer into a private namespace (to prevent conflicts)
55namespace {
56# include "AliFMDESDRevertexer.h"
57}
58
68aba90a 59
02a27b50 60class AliRawReader;
4347b38f 61
62//____________________________________________________________________
925e6570 63ClassImp(AliFMDReconstructor)
1a1fdef7 64#if 0
65 ; // This is here to keep Emacs for indenting the next line
66#endif
37c55dc0 67
4347b38f 68//____________________________________________________________________
69AliFMDReconstructor::AliFMDReconstructor()
70 : AliReconstructor(),
e9fd1e20 71 fMult(0x0),
72 fNMult(0),
73 fTreeR(0x0),
74 fCurrentVertex(0),
75 fESDObj(0x0),
9684be2f 76 fNoiseFactor(0),
77 fAngleCorrect(kTRUE),
78 fVertexType(kNoVertex),
9684be2f 79 fESD(0x0),
68aba90a 80 fDiagnostics(kTRUE),
9684be2f 81 fDiagStep1(0),
82 fDiagStep2(0),
83 fDiagStep3(0),
84 fDiagStep4(0),
85 fDiagAll(0)
4347b38f 86{
8f6ee336 87 // Make a new FMD reconstructor object - default CTOR.
a9579262 88 SetNoiseFactor();
89 SetAngleCorrect();
9b98d361 90 if (AliDebugLevel() > 0) fDiagnostics = kTRUE;
4347b38f 91}
92
42403906 93
94//____________________________________________________________________
95AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
8f6ee336 96 : AliReconstructor(),
97 fMult(other.fMult),
e9fd1e20 98 fNMult(other.fNMult),
99 fTreeR(other.fTreeR),
100 fCurrentVertex(other.fCurrentVertex),
101 fESDObj(other.fESDObj),
a9579262 102 fNoiseFactor(other.fNoiseFactor),
9684be2f 103 fAngleCorrect(other.fAngleCorrect),
104 fVertexType(other.fVertexType),
9684be2f 105 fESD(other.fESD),
106 fDiagnostics(other.fDiagnostics),
107 fDiagStep1(other.fDiagStep1),
108 fDiagStep2(other.fDiagStep2),
109 fDiagStep3(other.fDiagStep3),
110 fDiagStep4(other.fDiagStep4),
111 fDiagAll(other.fDiagAll)
42403906 112{
56b1929b 113 // Copy constructor
42403906 114}
115
116
117//____________________________________________________________________
118AliFMDReconstructor&
119AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
120{
56b1929b 121 // Assignment operator
a9579262 122 fMult = other.fMult;
123 fNMult = other.fNMult;
124 fTreeR = other.fTreeR;
e9fd1e20 125 fCurrentVertex = other.fCurrentVertex;
a9579262 126 fESDObj = other.fESDObj;
a9579262 127 fNoiseFactor = other.fNoiseFactor;
128 fAngleCorrect = other.fAngleCorrect;
9684be2f 129 fVertexType = other.fVertexType;
9684be2f 130 fESD = other.fESD;
131 fDiagnostics = other.fDiagnostics;
132 fDiagStep1 = other.fDiagStep1;
133 fDiagStep2 = other.fDiagStep2;
134 fDiagStep3 = other.fDiagStep3;
135 fDiagStep4 = other.fDiagStep4;
136 fDiagAll = other.fDiagAll;
42403906 137 return *this;
138}
56b1929b 139
140//____________________________________________________________________
141AliFMDReconstructor::~AliFMDReconstructor()
142{
143 // Destructor
8f6ee336 144 if (fMult) fMult->Delete();
145 if (fMult) delete fMult;
146 if (fESDObj) delete fESDObj;
56b1929b 147}
4347b38f 148
149//____________________________________________________________________
150void
d76c31f4 151AliFMDReconstructor::Init()
1a1fdef7 152{
153 // Initialize the reconstructor
f011bd38 154
bf000c32 155 // Initialize the geometry
9b48326f 156 AliFMDGeometry* geom = AliFMDGeometry::Instance();
157 geom->Init();
158 geom->InitTransformations();
159
160 // Initialize the parameters
161 AliFMDParameters* param = AliFMDParameters::Instance();
162 param->Init();
bf000c32 163
8f6ee336 164 // Current vertex position
1a1fdef7 165 fCurrentVertex = 0;
8f6ee336 166 // Create array of reconstructed strip multiplicities
bf000c32 167 fMult = new TClonesArray("AliFMDRecPoint", 51200);
8f6ee336 168 // Create ESD output object
169 fESDObj = new AliESDFMD;
170
9684be2f 171 // Check if we need diagnostics histograms
172 if (!fDiagnostics) return;
9b98d361 173 AliInfo("Making diagnostics histograms");
9684be2f 174 fDiagStep1 = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
175 1024, -.5, 1023.5, 1024, -.5, 1023.5);
176 fDiagStep1->SetDirectory(0);
177 fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
178 fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)",
179 fNoiseFactor));
180 fDiagStep2 = new TH2F("diagStep2", "ADC vs Edep deduced",
181 1024, -.5, 1023.5, 100, 0, 2);
182 fDiagStep2->SetDirectory(0);
183 fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
184 fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
185 fDiagStep3 = new TH2F("diagStep3", "Edep vs Edep path corrected",
186 100, 0., 2., 100, 0., 2.);
187 fDiagStep3->SetDirectory(0);
188 fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
189 fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
190 fDiagStep4 = new TH2F("diagStep4", "Edep vs Multiplicity deduced",
191 100, 0., 2., 100, -.1, 19.9);
192 fDiagStep4->SetDirectory(0);
193 fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
194 fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
195 fDiagAll = new TH2F("diagAll", "Read ADC vs Multiplicity deduced",
196 1024, -.5, 1023.5, 100, -.1, 19.9);
197 fDiagAll->SetDirectory(0);
198 fDiagAll->GetXaxis()->SetTitle("ADC (read)");
199 fDiagAll->GetYaxis()->SetTitle("Multiplicity");
4347b38f 200}
dc8af42e 201
4347b38f 202//____________________________________________________________________
203void
1a1fdef7 204AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
205 TTree* digitsTree) const
206{
207 // Convert Raw digits to AliFMDDigit's in a tree
f95a63c4 208 AliFMDDebug(2, ("Reading raw data into digits tree"));
1a1fdef7 209 AliFMDRawReader rawRead(reader, digitsTree);
210 // rawRead.SetSampleRate(fFMD->GetSampleRate());
211 rawRead.Exec();
5cf05dbb 212 AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
213 for (size_t i = 1; i <= 3; i++) {
214 fZS[i] = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
215 fZSFactor[i] = rawRead.NoiseFactor(map->Detector2DDL(i));
216 }
4347b38f 217}
218
4347b38f 219//____________________________________________________________________
220void
9684be2f 221AliFMDReconstructor::GetVertex() const
4347b38f 222{
97e94238 223 // Return the vertex to use.
224 // This is obtained from the ESD object.
225 // If not found, a warning is issued.
9684be2f 226 fVertexType = kNoVertex;
227 fCurrentVertex = 0;
a3537838 228 if (fESD) {
50b9d194 229 const AliESDVertex* vertex = fESD->GetPrimaryVertex();
f70f588a 230 if (!vertex) vertex = fESD->GetPrimaryVertexSPD();
231 if (!vertex) vertex = fESD->GetPrimaryVertexTPC();
232 if (!vertex) vertex = fESD->GetVertex();
50b9d194 233
8f6ee336 234 if (vertex) {
50b9d194 235 AliFMDDebug(2, ("Got %s (%s) from ESD: %f",
236 vertex->GetName(), vertex->GetTitle(), vertex->GetZv()));
8f6ee336 237 fCurrentVertex = vertex->GetZv();
9684be2f 238 fVertexType = kESDVertex;
239 return;
8f6ee336 240 }
50b9d194 241 else if (fESD->GetESDTZERO()) {
242 AliFMDDebug(2, ("Got primary vertex from T0: %f", fESD->GetT0zVertex()));
243 fCurrentVertex = fESD->GetT0zVertex();
244 fVertexType = kESDVertex;
245 return;
246 }
a3537838 247 }
9684be2f 248 AliWarning("Didn't get any vertex from ESD or generator");
249}
250
251
252//____________________________________________________________________
253void
50b9d194 254AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
ddaa8027 255{
256 // Reconstruct directly from raw data (no intermediate output on
257 // digit tree or rec point tree).
50b9d194 258 //
ddaa8027 259 // Parameters:
260 // reader Raw event reader
261 // ctree Not used.
50b9d194 262 AliFMDRawReader rawReader(reader, 0);
263
264 UShort_t det, sec, str, fac;
265 Short_t adc, oldDet = -1;
266 Bool_t zs;
267 Char_t rng;
268
269 while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
270 if (det != oldDet) {
271 fZS[det-1] = zs;
272 fZSFactor[det-1] = fac;
273 oldDet = det;
274 }
275 ProcessSignal(det, rng, sec, str, adc);
276 }
ddaa8027 277}
278
279//____________________________________________________________________
280void
9684be2f 281AliFMDReconstructor::Reconstruct(TTree* digitsTree,
282 TTree* clusterTree) const
283{
284 // Reconstruct event from digits in tree
285 // Get the FMD branch holding the digits.
286 // FIXME: The vertex may not be known yet, so we may have to move
287 // some of this to FillESD.
50b9d194 288 //
289 // Parameters:
290 // digitsTree Pointer to a tree containing digits
291 // clusterTree Pointer to output tree
292 //
f95a63c4 293 AliFMDDebug(2, ("Reconstructing from digits in a tree"));
9684be2f 294 GetVertex();
295
1a1fdef7 296 TBranch *digitBranch = digitsTree->GetBranch("FMD");
4347b38f 297 if (!digitBranch) {
1a1fdef7 298 Error("Exec", "No digit branch for the FMD found");
299 return;
4347b38f 300 }
0abe7182 301 TClonesArray* digits = new TClonesArray("AliFMDDigit");
1a1fdef7 302 digitBranch->SetAddress(&digits);
4347b38f 303
8f6ee336 304 if (fMult) fMult->Clear();
305 if (fESDObj) fESDObj->Clear();
306
307 fNMult = 0;
308 fTreeR = clusterTree;
309 fTreeR->Branch("FMD", &fMult);
310
f95a63c4 311 AliFMDDebug(5, ("Getting entry 0 from digit branch"));
1a1fdef7 312 digitBranch->GetEntry(0);
313
f95a63c4 314 AliFMDDebug(5, ("Processing digits"));
e802be3e 315 ProcessDigits(digits);
316
8f6ee336 317 Int_t written = clusterTree->Fill();
f95a63c4 318 AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
0abe7182 319 digits->Delete();
320 delete digits;
4347b38f 321}
1a1fdef7 322
8f6ee336 323
4347b38f 324//____________________________________________________________________
e802be3e 325void
326AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
4347b38f 327{
69b696b9 328 // For each digit, find the pseudo rapdity, azimuthal angle, and
329 // number of corrected ADC counts, and pass it on to the algorithms
330 // used.
50b9d194 331 //
332 // Parameters:
333 // digits Array of digits
334 //
e802be3e 335 Int_t nDigits = digits->GetEntries();
f95a63c4 336 AliFMDDebug(1, ("Got %d digits", nDigits));
a9579262 337 fESDObj->SetNoiseFactor(fNoiseFactor);
338 fESDObj->SetAngleCorrected(fAngleCorrect);
e802be3e 339 for (Int_t i = 0; i < nDigits; i++) {
340 AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
50b9d194 341 if (!digit) continue;
342 ProcessDigit(digit);
343 }
344}
8f6ee336 345
50b9d194 346//____________________________________________________________________
347void
348AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
349{
350 UShort_t det = digit->Detector();
351 Char_t rng = digit->Ring();
352 UShort_t sec = digit->Sector();
353 UShort_t str = digit->Strip();
354 Short_t adc = digit->Counts();
355
356 ProcessSignal(det, rng, sec, str, adc);
357}
358
359//____________________________________________________________________
360void
361AliFMDReconstructor::ProcessSignal(UShort_t det,
362 Char_t rng,
363 UShort_t sec,
364 UShort_t str,
365 Short_t adc) const
366{
367 // Process the signal from a single strip
368 //
369 // Parameters:
370 // det Detector ID
371 // rng Ring ID
372 // sec Sector ID
373 // rng Strip ID
374 // adc ADC counts
375 //
376 AliFMDParameters* param = AliFMDParameters::Instance();
377 // Check that the strip is not marked as dead
378 if (param->IsDead(det, rng, sec, str)) {
379 AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
380 return;
381 }
382
383 // digit->Print();
384 // Get eta and phi
385 Float_t eta, phi;
386 PhysicalCoordinates(det, rng, sec, str, eta, phi);
4347b38f 387
50b9d194 388 // Substract pedestal.
389 UShort_t counts = SubtractPedestal(det, rng, sec, str, adc);
390 if(counts == USHRT_MAX) return;
391
8f6ee336 392 // Gain match digits.
50b9d194 393 Double_t edep = Adc2Energy(det, rng, sec, str, eta, counts);
394 // Get rid of nonsense energy
395 if(edep < 0) return;
396
397 // Make rough multiplicity
398 Double_t mult = Energy2Multiplicity(det, rng, sec, str, edep);
399 // Get rid of nonsense mult
400 if (mult < 0) return;
401 AliFMDDebug(5, ("FMD%d%c[%2d,%3d]: "
402 "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
403 det, rng, sec, str, adc, counts, edep, mult));
404
405 // Create a `RecPoint' on the output branch.
406 if (fMult) {
407 AliFMDRecPoint* m =
408 new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
409 eta, phi, edep, mult);
410 (void)m; // Suppress warnings about unused variables.
411 fNMult++;
4347b38f 412 }
50b9d194 413
414 fESDObj->SetMultiplicity(det, rng, sec, str, mult);
415 fESDObj->SetEta(det, rng, sec, str, eta);
416
417 if (fDiagAll) fDiagAll->Fill(adc, mult);
418
4347b38f 419}
8f6ee336 420
50b9d194 421
422
1a1fdef7 423//____________________________________________________________________
424UShort_t
50b9d194 425AliFMDReconstructor::SubtractPedestal(UShort_t det,
426 Char_t rng,
427 UShort_t sec,
428 UShort_t str,
429 Short_t adc) const
1a1fdef7 430{
431 // Member function to subtract the pedestal from a digit
50b9d194 432 //
433 // Parameters:
434 // det Detector ID
435 // rng Ring ID
436 // sec Sector ID
437 // rng Strip ID
438 // adc # of ADC counts
439 // Return:
440 // Pedestal subtracted signal or USHRT_MAX in case of problems
441 //
8f6ee336 442 AliFMDParameters* param = AliFMDParameters::Instance();
50b9d194 443 Bool_t zs = fZS[det-1];
444 UShort_t fac = fZSFactor[det-1];
5cf05dbb 445 Float_t ped = (zs ? 0 :
50b9d194 446 param->GetPedestal(det, rng, sec, str));
447 Float_t noise = param->GetPedestalWidth(det, rng, sec, str);
68aba90a 448 if(ped < 0 || noise < 0) {
50b9d194 449 AliWarning(Form("Invalid pedestal (%f) or noise (%f) "
450 "for FMD%d%c[%02d,%03d]", ped, noise, det, rng, sec, str));
68aba90a 451 return USHRT_MAX;
452 }
453
50b9d194 454 AliFMDDebug(5, ("Subtracting pedestal %f from signal %d", ped, adc));
2aeec17d 455 // if (digit->Count3() > 0) adc = digit->Count3();
456 // else if (digit->Count2() > 0) adc = digit->Count2();
457 // else adc = digit->Count1();
5cf05dbb 458 Int_t counts = adc + Int_t(zs ? fac * noise : - ped);
459 counts = TMath::Max(Int_t(counts), 0);
a9579262 460 if (counts < noise * fNoiseFactor) counts = 0;
f95a63c4 461 if (counts > 0) AliFMDDebug(15, ("Got a hit strip"));
9684be2f 462 if (fDiagStep1) fDiagStep1->Fill(adc, counts);
8f6ee336 463
1a1fdef7 464 return UShort_t(counts);
465}
466
4347b38f 467//____________________________________________________________________
8f6ee336 468Float_t
50b9d194 469AliFMDReconstructor::Adc2Energy(UShort_t det,
470 Char_t rng,
471 UShort_t sec,
472 UShort_t str,
473 Float_t eta,
474 UShort_t count) const
8f6ee336 475{
476 // Converts number of ADC counts to energy deposited.
477 // Note, that this member function can be overloaded by derived
478 // classes to do strip-specific look-ups in databases or the like,
479 // to find the proper gain for a strip.
480 //
68aba90a 481 // In the first simple version, we calculate the energy deposited as
8f6ee336 482 //
483 // EnergyDeposited = cos(theta) * gain * count
484 //
485 // where
486 //
487 // Pre_amp_MIP_Range
488 // gain = ----------------- * Energy_deposited_per_MIP
489 // ADC_channel_size
490 //
491 // is constant and the same for all strips.
68aba90a 492 //
493 // For the production we use the conversion measured in the NBI lab.
494 // The total conversion is then:
50b9d194 495 //
496 // gain = ADC / DAC
497 //
498 // EdepMip * count
499 // => energy = ----------------
500 // gain * DACPerADC
501 //
502 // Parameters:
503 // det Detector ID
504 // rng Ring ID
505 // sec Sector ID
506 // rng Strip ID
507 // eta Psuedo-rapidity
508 // counts Number of ADC counts over pedestal
509 // Return
510 // The energy deposited in a single strip, or -1 in case of problems
68aba90a 511 //
a9579262 512 if (count <= 0) return 0;
8f6ee336 513 AliFMDParameters* param = AliFMDParameters::Instance();
50b9d194 514 Float_t gain = param->GetPulseGain(det, rng, sec, str);
68aba90a 515 // 'Tagging' bad gains as bad energy
516 if (gain < 0) {
50b9d194 517 AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
518 gain, det, rng, sec, str));
68aba90a 519 return -1;
520 }
50b9d194 521 AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
68aba90a 522 count, gain,param->GetDACPerMIP()));
a9579262 523
50b9d194 524 Double_t edep = ((count * param->GetEdepMip())
525 / (gain * param->GetDACPerMIP()));
9684be2f 526 if (fDiagStep2) fDiagStep2->Fill(count, edep);
a9579262 527 if (fAngleCorrect) {
9684be2f 528 Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
529 Double_t corr = TMath::Abs(TMath::Cos(theta));
530 Double_t cedep = corr * edep;
f95a63c4 531 AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
9684be2f 532 edep, corr, cedep, eta, theta));
533 if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
534 edep = cedep;
a9579262 535 }
8f6ee336 536 return edep;
537}
538
539//____________________________________________________________________
540Float_t
50b9d194 541AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
542 Char_t /*rng*/,
543 UShort_t /*sec*/,
544 UShort_t /*str*/,
545 Float_t edep) const
8f6ee336 546{
547 // Converts an energy signal to number of particles.
548 // Note, that this member function can be overloaded by derived
549 // classes to do strip-specific look-ups in databases or the like,
550 // to find the proper gain for a strip.
551 //
552 // In this simple version, we calculate the multiplicity as
553 //
554 // multiplicity = Energy_deposited / Energy_deposited_per_MIP
555 //
556 // where
557 //
558 // Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
559 //
560 // is constant and the same for all strips
50b9d194 561 //
562 // Parameters:
563 // det Detector ID
564 // rng Ring ID
565 // sec Sector ID
566 // rng Strip ID
567 // edep Energy deposited in a single strip
568 // Return
569 // The "bare" multiplicity corresponding to the energy deposited
8f6ee336 570 AliFMDParameters* param = AliFMDParameters::Instance();
571 Double_t edepMIP = param->GetEdepMip();
572 Float_t mult = edep / edepMIP;
68aba90a 573#if 0
8f6ee336 574 if (edep > 0)
f95a63c4 575 AliFMDDebug(15, ("Translating energy %f to multiplicity via "
8f6ee336 576 "divider %f->%f", edep, edepMIP, mult));
68aba90a 577#endif
9684be2f 578 if (fDiagStep4) fDiagStep4->Fill(edep, mult);
8f6ee336 579 return mult;
580}
581
582//____________________________________________________________________
583void
50b9d194 584AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
585 Char_t rng,
586 UShort_t sec,
587 UShort_t str,
8f6ee336 588 Float_t& eta,
589 Float_t& phi) const
590{
591 // Get the eta and phi of a digit
592 //
50b9d194 593 // Parameters:
594 // det Detector ID
595 // rng Ring ID
596 // sec Sector ID
597 // rng Strip ID
598 // eta On return, contains the psuedo-rapidity of the strip
599 // phi On return, contains the azimuthal angle of the strip
600 //
9b48326f 601 AliFMDGeometry* geom = AliFMDGeometry::Instance();
bf000c32 602 Double_t x, y, z, r, theta;
50b9d194 603 geom->Detector2XYZ(det, rng, sec, str, x, y, z);
bf000c32 604 // Correct for vertex offset.
605 z += fCurrentVertex;
606 phi = TMath::ATan2(y, x);
607 r = TMath::Sqrt(y * y + x * x);
608 theta = TMath::ATan2(r, z);
609 eta = -TMath::Log(TMath::Tan(theta / 2));
8f6ee336 610}
611
612
613
614//____________________________________________________________________
4347b38f 615void
1a1fdef7 616AliFMDReconstructor::FillESD(TTree* /* digitsTree */,
617 TTree* /* clusterTree */,
af885e0f 618 AliESDEvent* esd) const
121a60bd 619{
42403906 620 // nothing to be done
69b696b9 621 // FIXME: The vertex may not be known when Reconstruct is executed,
622 // so we may have to move some of that member function here.
f95a63c4 623 AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
8f6ee336 624 // fESDObj->Print();
625
50b9d194 626 Double_t oldVz = fCurrentVertex;
627 GetVertex();
628 if (fVertexType != kNoVertex) {
629 AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
630 fCurrentVertex, oldVz));
631 AliFMDESDRevertexer revertexer;
632 revertexer.Revertex(fESDObj, fCurrentVertex);
633 }
634
8f6ee336 635 if (esd) {
f95a63c4 636 AliFMDDebug(2, ("Writing FMD data to ESD tree"));
8f6ee336 637 esd->SetFMDData(fESDObj);
a3537838 638 }
9684be2f 639
640 if (!fDiagnostics || !esd) return;
641 static bool first = true;
69893a66 642 // This is most likely NOT the event number you'd like to use. It
643 // has nothing to do with the 'real' event number.
644 // - That's OK. We just use it for the name of the directory -
645 // nothing else. Christian
646 Int_t evno = esd->GetEventNumberInFile();
647 AliFMDDebug(1, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
9684be2f 648 TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
649 first = false;
650 f.cd();
651 TDirectory* d = f.mkdir(Form("%03d", evno),
652 Form("Diagnostics histograms for event # %d", evno));
653 d->cd();
654 if (fDiagStep1) fDiagStep1->Write();
655 if (fDiagStep2) fDiagStep2->Write();
656 if (fDiagStep3) fDiagStep3->Write();
657 if (fDiagStep4) fDiagStep4->Write();
658 if (fDiagAll) fDiagAll->Write();
659 d->Write();
660 f.Write();
661 f.Close();
662
663 if (fDiagStep1) fDiagStep1->Reset();
664 if (fDiagStep2) fDiagStep2->Reset();
665 if (fDiagStep3) fDiagStep3->Reset();
666 if (fDiagStep4) fDiagStep4->Reset();
667 if (fDiagAll) fDiagAll->Reset();
121a60bd 668}
669
42403906 670//____________________________________________________________________
ddaa8027 671void
672AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
673 AliESDEvent* esd) const
674{
675 TTree* dummy = 0;
676 FillESD(dummy, clusterTree, esd);
677}
678
679//____________________________________________________________________
42403906 680//
681// EOF
682//