Added docs and fixed a bug
[u/mrichter/AliRoot.git] / FMD / AliFMDMultPoisson.cxx
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 // 
20 // Reconstruct charged particle multiplicity in the FMD 
21 // 
22 // [See also the AliFMDReconstructor class]
23 // 
24 // This class reconstructs the muliplicity in regions based on the
25 // ratio of empty to full strips. 
26 //
27 #include "AliFMD.h"                     // ALIFMD_H
28 #include "AliFMDMultPoisson.h"          // ALIFMDMULTPOISSON_H
29 #include "AliFMDMultRegion.h"           // ALIFMDMULTREGION_H
30 #include "AliFMDDigit.h"                // ALIFMDDIGIT_H
31 #include "AliLog.h"                     // ALILOG_H
32 #include <TClonesArray.h>               // ROOT_TClonesArray
33 #include <TTree.h>                      // ROOT_TTree
34
35 //____________________________________________________________________
36 ClassImp(AliFMDMultPoisson);
37
38 //____________________________________________________________________
39 AliFMDMultPoisson::AliFMDMultPoisson()
40   : AliFMDMultAlgorithm("Poisson", "Poisson"),
41     fDeltaEta(0), 
42     fDeltaPhi(0), 
43     fThreshold(0)
44 {
45   SetDeltaEta();
46   SetDeltaPhi();
47   SetThreshold();
48   fMult = new TClonesArray("AliFMDMultRegion", 1000);
49 }
50
51 //____________________________________________________________________
52 void
53 AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ) 
54 {
55   // Reset internal data
56   AliFMDMultAlgorithm::PreEvent(tree, ipZ);
57   fCurrentVertexZ = ipZ;
58   fEmpty.Clear(kFALSE);
59
60   // Make a branch in the reconstruction tree. 
61   const Int_t kBufferSize = 16000;
62   fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);  
63   
64 }
65
66 //____________________________________________________________________
67 void
68 AliFMDMultPoisson::ProcessDigit(AliFMDDigit*  digit, 
69                                 Float_t       /* eta */, 
70                                 Float_t       /* phi */, 
71                                 UShort_t      count)
72 {
73   // Process one digit. 
74   // 
75   // Parameters: 
76   //    
77   //   digit            Digit to process 
78   //   ipZ              Z--coordinate of the primary interaction
79   //                    vertex of this event 
80   //
81   if (!digit) return;
82   if (count < fThreshold) fEmpty(digit->Detector() - 1, 
83                                  digit->Ring(), 
84                                  digit->Sector(), 
85                                  digit->Strip()) = kTRUE;
86 }
87
88 //____________________________________________________________________
89 void
90 AliFMDMultPoisson::PostEvent() 
91 {
92   // Fill the branch 
93   // Based on the information in the cache, do the reconstruction. 
94
95   // Loop over the detectors 
96   for (Int_t i = 1; i <= 3; i++) {
97     AliFMDSubDetector* sub = 0;
98     switch (i) {
99     case 1: sub = fFMD->GetFMD1(); break;
100     case 2: sub = fFMD->GetFMD2(); break;
101     case 3: sub = fFMD->GetFMD3(); break;
102     }
103     if (!sub) continue;
104         
105     // Loop over the rings in the detector
106     for (Int_t j = 0; j < 2; j++) {
107       Float_t     rZ = 0;
108       AliFMDRing* r  = 0;
109       switch (j) {
110       case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
111       case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
112       }
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         
149         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
150                           minPhi, maxPhi, minSector, maxSector));
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));
161           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
162           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
163
164           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
165                             "    [radii %f, %f, thetas %f, %f, sign %d]", 
166                             minEta, maxEta, minStrip, maxStrip, 
167                             minR, maxR, theta1, theta2, sign));
168
169           // Count number of empty strips
170           Int_t   emptyStrips = 0;
171           for (Int_t sector = minSector; sector < maxSector; sector++) 
172             for (Int_t strip = minStrip; strip < maxStrip; strip++) 
173               if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip)) 
174                 emptyStrips++;
175           
176           // The total number of strips 
177           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
178           
179           // Log ratio of empty to total number of strips 
180           AliDebug(10, Form("Lambda= %d / %d = %f", 
181                             emptyStrips, nTotal, 
182                             Float_t(emptyStrips) / nTotal));
183           
184           Double_t lambda = (emptyStrips > 0 ? 
185                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
186                              1);
187
188           // The reconstructed number of particles is then given by 
189           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
190             
191           // Add a AliFMDMultRegion to the reconstruction tree. 
192           AliFMDMultRegion* m = new((*fMult)[fNMult])   
193             AliFMDMultRegion(sub->GetId(), r->GetId(),
194                              minSector, maxSector, minStrip, maxStrip,
195                              minEta, maxEta, minPhi, maxPhi,
196                              reconstructed, AliFMDMultRegion::kPoission);
197           (void)m;
198           fNMult++;
199         } // phi 
200       } // eta
201     } // ring 
202   } // detector 
203 }
204
205
206 //____________________________________________________________________
207 // 
208 // EOF
209 //