Fixed some coding style violations.
[u/mrichter/AliRoot.git] / FMD / AliFMDMultPoisson.cxx
CommitLineData
56b1929b 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
16/* $Id$ */
17
18//____________________________________________________________________
19//
7684b53c 20// Reconstruct charged particle multiplicity in the FMD
21//
22// [See also the AliFMDReconstructor class]
56b1929b 23//
7684b53c 24// This class reconstructs the muliplicity in regions based on the
25// ratio of empty to full strips.
088f8e79 26// Hence the name `poisson'
7684b53c 27//
56b1929b 28#include "AliFMD.h" // ALIFMD_H
1a1fdef7 29#include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
30#include "AliFMDDetector.h" // ALIFMDDETECTOR_H
31#include "AliFMDRing.h" // ALIFMDRING_H
56b1929b 32#include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
33#include "AliFMDMultRegion.h" // ALIFMDMULTREGION_H
34#include "AliFMDDigit.h" // ALIFMDDIGIT_H
35#include "AliLog.h" // ALILOG_H
36#include <TClonesArray.h> // ROOT_TClonesArray
37#include <TTree.h> // ROOT_TTree
38
39//____________________________________________________________________
925e6570 40ClassImp(AliFMDMultPoisson)
1a1fdef7 41#if 0
42 ; // This is here to keep Emacs for indenting the next line
43#endif
56b1929b 44
45//____________________________________________________________________
46AliFMDMultPoisson::AliFMDMultPoisson()
47 : AliFMDMultAlgorithm("Poisson", "Poisson"),
48 fDeltaEta(0),
49 fDeltaPhi(0),
50 fThreshold(0)
51{
088f8e79 52 // CTOR
56b1929b 53 SetDeltaEta();
54 SetDeltaPhi();
55 SetThreshold();
56 fMult = new TClonesArray("AliFMDMultRegion", 1000);
57}
58
59//____________________________________________________________________
60void
61AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ)
62{
63 // Reset internal data
64 AliFMDMultAlgorithm::PreEvent(tree, ipZ);
65 fCurrentVertexZ = ipZ;
69b696b9 66 fEmpty.Reset(kFALSE);
56b1929b 67
68 // Make a branch in the reconstruction tree.
69 const Int_t kBufferSize = 16000;
70 fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);
71
72}
73
74//____________________________________________________________________
75void
76AliFMDMultPoisson::ProcessDigit(AliFMDDigit* digit,
77 Float_t /* eta */,
78 Float_t /* phi */,
79 UShort_t count)
80{
81 // Process one digit.
82 //
83 // Parameters:
84 //
85 // digit Digit to process
86 // ipZ Z--coordinate of the primary interaction
87 // vertex of this event
88 //
89 if (!digit) return;
90 if (count < fThreshold) fEmpty(digit->Detector() - 1,
91 digit->Ring(),
92 digit->Sector(),
93 digit->Strip()) = kTRUE;
94}
95
96//____________________________________________________________________
97void
98AliFMDMultPoisson::PostEvent()
99{
100 // Fill the branch
101 // Based on the information in the cache, do the reconstruction.
102
103 // Loop over the detectors
104 for (Int_t i = 1; i <= 3; i++) {
1a1fdef7 105 AliFMDGeometry* fmd = AliFMDGeometry::Instance();
106 AliFMDDetector* sub = fmd->GetDetector(i);
56b1929b 107 if (!sub) continue;
108
109 // Loop over the rings in the detector
110 for (Int_t j = 0; j < 2; j++) {
1a1fdef7 111 AliFMDRing* r = sub->GetRing((j == 0 ? 'I' : 'O'));
112 Float_t rZ = sub->GetRingZ((j == 0 ? 'I' : 'O'));
56b1929b 113 if (!r) continue;
114
115 // Calculate low/high theta and eta
116 // FIXME: Is this right?
117 Float_t realZ = fCurrentVertexZ + rZ;
118 Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
119 Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
120 Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
121 Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
122 if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
123 Float_t tmp = etaIn;
124 etaIn = etaOut;
125 etaOut = tmp;
126 }
127
128 //-------------------------------------------------------------
129 //
130 // Here starts poisson method
131 //
132 // Calculate eta step per strip, number of eta steps, number of
133 // phi steps, and check the sign of the eta increment
134 Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
135 Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
136 Int_t nPhi = Int_t(360. / fDeltaPhi);
137 Float_t sign = TMath::Sign(Float_t(1.), etaIn);
138
139 AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
140 sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
141
142 // Loop over relevant phi values
143 for (Int_t p = 0; p < nPhi; p++) {
144 Float_t minPhi = p * fDeltaPhi;
145 Float_t maxPhi = minPhi + fDeltaPhi;
146 UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
147 UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
148
c263f0ab 149 // AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
150 // minPhi, maxPhi, minSector, maxSector));
56b1929b 151 // Loop over relevant eta values
152 for (Int_t e = nEta; e >= 0; --e) {
153 Float_t maxEta = etaIn - sign * e * fDeltaEta;
154 Float_t minEta = maxEta - sign * fDeltaEta;
155 if (sign > 0) minEta = TMath::Max(minEta, etaOut);
156 else minEta = TMath::Min(minEta, etaOut);
157 Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
158 Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
159 Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
160 Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
a3537838 161 // Calculate the weighted mean eta of the region
162 Float_t minW2 = TMath::Power(minR * 2 * TMath::Pi() *
163 ((maxPhi - minPhi)/360),2);
164 Float_t maxW2 = TMath::Power(minR * 2 * TMath::Pi() *
165 ((maxPhi - minPhi)/360), 2);
166 Float_t meanEta = ((minEta / minW2 + maxEta / maxW2) /
167 (1 / (minW2 + maxW2)));
c263f0ab 168 //UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
169 // UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
170
a3537838 171 UShort_t minStrip = UShort_t(r->GetNStrips() -
172 (etaIn - minEta) * stripEta + 0.5);
173 UShort_t maxStrip = UShort_t(r->GetNStrips() -
174 (etaIn - maxEta) * stripEta + 0.5);
56b1929b 175
176 AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
c263f0ab 177 " [radii %f, %f, thetas %f, %f, sign %d]",
178 minEta, maxEta, minStrip, maxStrip,
56b1929b 179 minR, maxR, theta1, theta2, sign));
c263f0ab 180
56b1929b 181 // Count number of empty strips
182 Int_t emptyStrips = 0;
183 for (Int_t sector = minSector; sector < maxSector; sector++)
184 for (Int_t strip = minStrip; strip < maxStrip; strip++)
185 if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip))
186 emptyStrips++;
187
188 // The total number of strips
189 Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
56b1929b 190 // Log ratio of empty to total number of strips
56b1929b 191
192 Double_t lambda = (emptyStrips > 0 ?
193 - TMath::Log(Double_t(emptyStrips) / nTotal) :
194 1);
195
196 // The reconstructed number of particles is then given by
197 Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
c263f0ab 198 AliDebug(10, Form("Lambda= %d / %f = %f particles %d",
199 emptyStrips, nTotal,
200 Float_t(emptyStrips) / nTotal , reconstructed));
56b1929b 201
202 // Add a AliFMDMultRegion to the reconstruction tree.
7684b53c 203 AliFMDMultRegion* m = new((*fMult)[fNMult])
56b1929b 204 AliFMDMultRegion(sub->GetId(), r->GetId(),
205 minSector, maxSector, minStrip, maxStrip,
a3537838 206 minEta, maxEta, meanEta, minPhi, maxPhi,
56b1929b 207 reconstructed, AliFMDMultRegion::kPoission);
7684b53c 208 (void)m;
56b1929b 209 fNMult++;
210 } // phi
211 } // eta
212 } // ring
213 } // detector
214}
215
216
217//____________________________________________________________________
218//
219// EOF
220//