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