1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //____________________________________________________________________
20 // Reconstruct charged particle multiplicity in the FMD
22 // [See also the AliFMDReconstructor class]
24 // This class reconstructs the muliplicity in regions based on the
25 // ratio of empty to full strips.
27 #include "AliFMD.h" // ALIFMD_H
28 #include "AliFMDGeometry.h" // ALIFMDGEOMETRY_H
29 #include "AliFMDDetector.h" // ALIFMDDETECTOR_H
30 #include "AliFMDRing.h" // ALIFMDRING_H
31 #include "AliFMDMultPoisson.h" // ALIFMDMULTPOISSON_H
32 #include "AliFMDMultRegion.h" // ALIFMDMULTREGION_H
33 #include "AliFMDDigit.h" // ALIFMDDIGIT_H
34 #include "AliLog.h" // ALILOG_H
35 #include <TClonesArray.h> // ROOT_TClonesArray
36 #include <TTree.h> // ROOT_TTree
38 //____________________________________________________________________
39 ClassImp(AliFMDMultPoisson)
41 ; // This is here to keep Emacs for indenting the next line
44 //____________________________________________________________________
45 AliFMDMultPoisson::AliFMDMultPoisson()
46 : AliFMDMultAlgorithm("Poisson", "Poisson"),
54 fMult = new TClonesArray("AliFMDMultRegion", 1000);
57 //____________________________________________________________________
59 AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ)
61 // Reset internal data
62 AliFMDMultAlgorithm::PreEvent(tree, ipZ);
63 fCurrentVertexZ = ipZ;
66 // Make a branch in the reconstruction tree.
67 const Int_t kBufferSize = 16000;
68 fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);
72 //____________________________________________________________________
74 AliFMDMultPoisson::ProcessDigit(AliFMDDigit* digit,
83 // digit Digit to process
84 // ipZ Z--coordinate of the primary interaction
85 // vertex of this event
88 if (count < fThreshold) fEmpty(digit->Detector() - 1,
91 digit->Strip()) = kTRUE;
94 //____________________________________________________________________
96 AliFMDMultPoisson::PostEvent()
99 // Based on the information in the cache, do the reconstruction.
101 // Loop over the detectors
102 for (Int_t i = 1; i <= 3; i++) {
103 AliFMDGeometry* fmd = AliFMDGeometry::Instance();
104 AliFMDDetector* sub = fmd->GetDetector(i);
107 // Loop over the rings in the detector
108 for (Int_t j = 0; j < 2; j++) {
109 AliFMDRing* r = sub->GetRing((j == 0 ? 'I' : 'O'));
110 Float_t rZ = sub->GetRingZ((j == 0 ? 'I' : 'O'));
113 // Calculate low/high theta and eta
114 // FIXME: Is this right?
115 Float_t realZ = fCurrentVertexZ + rZ;
116 Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
117 Float_t thetaIn = TMath::ATan2(r->GetLowR(), realZ);
118 Float_t etaOut = - TMath::Log(TMath::Tan(thetaOut / 2));
119 Float_t etaIn = - TMath::Log(TMath::Tan(thetaIn / 2));
120 if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
126 //-------------------------------------------------------------
128 // Here starts poisson method
130 // Calculate eta step per strip, number of eta steps, number of
131 // phi steps, and check the sign of the eta increment
132 Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
133 Int_t nEta = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
134 Int_t nPhi = Int_t(360. / fDeltaPhi);
135 Float_t sign = TMath::Sign(Float_t(1.), etaIn);
137 AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
138 sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
140 // Loop over relevant phi values
141 for (Int_t p = 0; p < nPhi; p++) {
142 Float_t minPhi = p * fDeltaPhi;
143 Float_t maxPhi = minPhi + fDeltaPhi;
144 UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
145 UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
147 AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
148 minPhi, maxPhi, minSector, maxSector));
149 // Loop over relevant eta values
150 for (Int_t e = nEta; e >= 0; --e) {
151 Float_t maxEta = etaIn - sign * e * fDeltaEta;
152 Float_t minEta = maxEta - sign * fDeltaEta;
153 if (sign > 0) minEta = TMath::Max(minEta, etaOut);
154 else minEta = TMath::Min(minEta, etaOut);
155 Float_t theta1 = 2 * TMath::ATan(TMath::Exp(-minEta));
156 Float_t theta2 = 2 * TMath::ATan(TMath::Exp(-maxEta));
157 Float_t minR = TMath::Abs(realZ * TMath::Tan(theta2));
158 Float_t maxR = TMath::Abs(realZ * TMath::Tan(theta1));
159 UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
160 UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
162 AliDebug(10, Form(" Now in eta range %f, %f (strips %d, %d)\n"
163 " [radii %f, %f, thetas %f, %f, sign %d]",
164 minEta, maxEta, minStrip, maxStrip,
165 minR, maxR, theta1, theta2, sign));
167 // Count number of empty strips
168 Int_t emptyStrips = 0;
169 for (Int_t sector = minSector; sector < maxSector; sector++)
170 for (Int_t strip = minStrip; strip < maxStrip; strip++)
171 if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip))
174 // The total number of strips
175 Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
177 // Log ratio of empty to total number of strips
178 AliDebug(10, Form("Lambda= %d / %d = %f",
180 Float_t(emptyStrips) / nTotal));
182 Double_t lambda = (emptyStrips > 0 ?
183 - TMath::Log(Double_t(emptyStrips) / nTotal) :
186 // The reconstructed number of particles is then given by
187 Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
189 // Add a AliFMDMultRegion to the reconstruction tree.
190 AliFMDMultRegion* m = new((*fMult)[fNMult])
191 AliFMDMultRegion(sub->GetId(), r->GetId(),
192 minSector, maxSector, minStrip, maxStrip,
193 minEta, maxEta, minPhi, maxPhi,
194 reconstructed, AliFMDMultRegion::kPoission);
204 //____________________________________________________________________