Fixes, and extra debug
[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 "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
37
38 //____________________________________________________________________
39 ClassImp(AliFMDMultPoisson)
40 #if 0
41   ; // This is here to keep Emacs for indenting the next line
42 #endif
43
44 //____________________________________________________________________
45 AliFMDMultPoisson::AliFMDMultPoisson()
46   : AliFMDMultAlgorithm("Poisson", "Poisson"),
47     fDeltaEta(0), 
48     fDeltaPhi(0), 
49     fThreshold(0)
50 {
51   SetDeltaEta();
52   SetDeltaPhi();
53   SetThreshold();
54   fMult = new TClonesArray("AliFMDMultRegion", 1000);
55 }
56
57 //____________________________________________________________________
58 void
59 AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ) 
60 {
61   // Reset internal data
62   AliFMDMultAlgorithm::PreEvent(tree, ipZ);
63   fCurrentVertexZ = ipZ;
64   fEmpty.Reset(kFALSE);
65
66   // Make a branch in the reconstruction tree. 
67   const Int_t kBufferSize = 16000;
68   fTreeR->Branch("FMDPoisson", &fMult, kBufferSize);  
69   
70 }
71
72 //____________________________________________________________________
73 void
74 AliFMDMultPoisson::ProcessDigit(AliFMDDigit*  digit, 
75                                 Float_t       /* eta */, 
76                                 Float_t       /* phi */, 
77                                 UShort_t      count)
78 {
79   // Process one digit. 
80   // 
81   // Parameters: 
82   //    
83   //   digit            Digit to process 
84   //   ipZ              Z--coordinate of the primary interaction
85   //                    vertex of this event 
86   //
87   if (!digit) return;
88   if (count < fThreshold) fEmpty(digit->Detector() - 1, 
89                                  digit->Ring(), 
90                                  digit->Sector(), 
91                                  digit->Strip()) = kTRUE;
92 }
93
94 //____________________________________________________________________
95 void
96 AliFMDMultPoisson::PostEvent() 
97 {
98   // Fill the branch 
99   // Based on the information in the cache, do the reconstruction. 
100
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);
105     if (!sub) continue;
106         
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'));
111       if (!r) continue;
112       
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)) {
121         Float_t tmp = etaIn;
122         etaIn       = etaOut;
123         etaOut      = tmp;
124       }
125
126       //-------------------------------------------------------------
127       //
128       // Here starts poisson method 
129       //
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);
136
137       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
138                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
139
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();
146         
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           // Calculate the weighted mean eta of the region
160           Float_t  minW2    = TMath::Power(minR * 2 * TMath::Pi() * 
161                                            ((maxPhi - minPhi)/360),2);
162           Float_t  maxW2    = TMath::Power(minR * 2 * TMath::Pi() * 
163                                            ((maxPhi - minPhi)/360), 2);
164           Float_t  meanEta  = ((minEta / minW2 + maxEta / maxW2) / 
165                                (1 / (minW2 + maxW2)));
166           //UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
167           // UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
168
169           UShort_t minStrip = UShort_t(r->GetNStrips() - 
170                                        (etaIn - minEta) * stripEta + 0.5);
171           UShort_t maxStrip = UShort_t(r->GetNStrips() - 
172                                        (etaIn - maxEta) * stripEta + 0.5);
173
174           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
175                             "    [radii %f, %f, thetas %f, %f, sign %d]", 
176                             minEta, maxEta, minStrip, maxStrip,
177                             minR, maxR, theta1, theta2, sign));
178           
179           // Count number of empty strips
180           Int_t   emptyStrips = 0;
181           for (Int_t sector = minSector; sector < maxSector; sector++) 
182             for (Int_t strip = minStrip; strip < maxStrip; strip++) 
183               if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip)) 
184                 emptyStrips++;
185           
186           // The total number of strips 
187           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
188           // Log ratio of empty to total number of strips 
189           
190           Double_t lambda = (emptyStrips > 0 ? 
191                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
192                              1);
193
194           // The reconstructed number of particles is then given by 
195           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
196           AliDebug(10, Form("Lambda= %d / %f = %f particles %d", 
197                             emptyStrips, nTotal, 
198                             Float_t(emptyStrips) / nTotal , reconstructed));
199             
200           // Add a AliFMDMultRegion to the reconstruction tree. 
201           AliFMDMultRegion* m = new((*fMult)[fNMult])   
202             AliFMDMultRegion(sub->GetId(), r->GetId(),
203                              minSector, maxSector, minStrip, maxStrip,
204                              minEta, maxEta, meanEta, minPhi, maxPhi,
205                              reconstructed, AliFMDMultRegion::kPoission);
206           (void)m;
207           fNMult++;
208         } // phi 
209       } // eta
210     } // ring 
211   } // detector 
212 }
213
214
215 //____________________________________________________________________
216 // 
217 // EOF
218 //