Line data Source code
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 : /* $Id$ */
16 : /** @file AliFMDv1.cxx
17 : @author Christian Holm Christensen <cholm@nbi.dk>
18 : @date Mon Mar 27 12:48:51 2006
19 : @brief Concrete implementation of FMD detector driver - detailed
20 : version
21 : @ingroup FMD_sim
22 : */
23 : //____________________________________________________________________
24 : //
25 : // Forward Multiplicity Detector based on Silicon wafers. This class
26 : // contains the base procedures for the Forward Multiplicity detector
27 : // Detector consists of 3 sub-detectors FMD1, FMD2, and FMD3, each of
28 : // which has 1 or 2 rings of silicon sensors.
29 : // This class contains the detailed version of the FMD - that is, hits
30 : // are produced during simulation.
31 : //
32 : // See also the class AliFMD for a more detailed explanation of the
33 : // various componets.
34 : //
35 : #include <TVirtualMC.h> // ROOT_TVirtualMC
36 : #include <AliRun.h> // ALIRUN_H
37 : #include <AliMC.h> // ALIMC_H
38 : // #include <AliLog.h> // ALILOG_H
39 : #include "AliFMDDebug.h" // Better debug macros
40 : #include "AliFMDv1.h" // ALIFMDV1_H
41 : // #include "AliFMDGeometryBuilder.h"
42 : #include "AliFMDGeometry.h"
43 : #include "AliFMDDetector.h"
44 : #include "AliFMDRing.h"
45 : #include <TParticlePDG.h>
46 : #include <TDatabasePDG.h>
47 : #include "AliFMDHit.h"
48 :
49 : //____________________________________________________________________
50 12 : ClassImp(AliFMDv1)
51 : #if 0
52 : ; // This is here to keep Emacs for indenting the next line
53 : #endif
54 :
55 :
56 : //____________________________________________________________________
57 : Bool_t
58 : AliFMDv1::VMC2FMD(TLorentzVector& v, UShort_t& detector,
59 : Char_t& ring, UShort_t& sector, UShort_t& strip) const
60 : {
61 : // Convert VMC coordinates to detector coordinates
62 0 : TVirtualMC* mc = TVirtualMC::GetMC();
63 0 : AliFMDGeometry* fmd = AliFMDGeometry::Instance();
64 :
65 : // Get track position
66 0 : mc->TrackPosition(v);
67 0 : Int_t moduleno; mc->CurrentVolOffID(fmd->GetModuleOff(), moduleno);
68 0 : Int_t iring; mc->CurrentVolOffID(fmd->GetRingOff(), iring);
69 0 : ring = Char_t(iring);
70 0 : Int_t det; mc->CurrentVolOffID(fmd->GetDetectorOff(), det);
71 0 : detector = det;
72 :
73 :
74 : // Get the ring geometry
75 : //Int_t nsec = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
76 0 : Int_t nstr = fmd->GetDetector(detector)->GetRing(ring)->GetNStrips();
77 0 : Double_t lowr = fmd->GetDetector(detector)->GetRing(ring)->GetMinR();
78 0 : Double_t theta = fmd->GetDetector(detector)->GetRing(ring)->GetTheta();
79 0 : Double_t pitch = fmd->GetDetector(detector)->GetRing(ring)->GetPitch();
80 :
81 : // Figure out the strip number
82 0 : Double_t r = TMath::Sqrt(v.X() * v.X() + v.Y() * v.Y());
83 0 : Int_t str = Int_t((r - lowr) / pitch);
84 0 : if (str < 0 || str >= nstr) return kFALSE;
85 0 : strip = str;
86 :
87 : // Figure out the sector number
88 0 : Double_t phi = TMath::ATan2(v.Y(), v.X()) * 180. / TMath::Pi();
89 0 : if (phi < 0) phi = 360. + phi;
90 0 : Double_t t = phi - 2 * moduleno * theta;
91 0 : sector = 2 * moduleno;
92 0 : if (t < 0 || t > 2 * theta) return kFALSE;
93 0 : else if (t > theta) sector += 1;
94 :
95 0 : AliFMDDebug(40, ("<1> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
96 : detector, ring, sector, strip, mc->CurrentVolPath()));
97 0 : return kTRUE;
98 0 : }
99 :
100 : //____________________________________________________________________
101 : Bool_t
102 : AliFMDv1::VMC2FMD(Int_t copy, TLorentzVector& v,
103 : UShort_t& detector, Char_t& ring,
104 : UShort_t& sector, UShort_t& strip) const
105 : {
106 : // Convert VMC coordinates to detector coordinates
107 676 : TVirtualMC* mc = TVirtualMC::GetMC();
108 338 : AliFMDGeometry* fmd = AliFMDGeometry::Instance();
109 :
110 338 : strip = copy - 1;
111 338 : Int_t sectordiv; mc->CurrentVolOffID(fmd->GetSectorOff(), sectordiv);
112 338 : if (fmd->GetModuleOff() >= 0) {
113 338 : Int_t module; mc->CurrentVolOffID(fmd->GetModuleOff(), module);
114 338 : sector = 2 * module + sectordiv;
115 338 : }
116 : else
117 0 : sector = sectordiv;
118 676 : AliFMDDebug(30, ("Getting ring volume with offset %d -> %s",
119 : fmd->GetRingOff(),
120 : mc->CurrentVolOffName(fmd->GetRingOff())));
121 338 : Int_t iring; mc->CurrentVolOffID(fmd->GetRingOff(), iring);
122 338 : ring = Char_t(iring);
123 338 : Int_t det; mc->CurrentVolOffID(fmd->GetDetectorOff(), det);
124 338 : detector = det;
125 :
126 : //Double_t rz = fmd->GetDetector(detector)->GetRingZ(ring);
127 338 : AliFMDDetector* gdet = fmd->GetDetector(detector);
128 338 : AliFMDRing* gring = gdet->GetRing(ring);
129 338 : if (!gring) {
130 0 : AliFatal(Form("Ring %c not found (volume was %s at offset %d in path %s)",
131 : ring,
132 : mc->CurrentVolOffName(fmd->GetRingOff()),
133 : fmd->GetRingOff(),
134 : mc->CurrentVolPath()));
135 0 : }
136 338 : Int_t n = gring->GetNSectors();
137 : #if 0
138 : if (rz < 0) {
139 : Int_t s = ((n - sector + n / 2) % n) + 1;
140 : AliFMDDebug(1, ("Recalculating sector to %d (=%d-%d+%d/2%%%d+1 z=%f)",
141 : s, n, sector, n, n, rz));
142 : sector = s;
143 : }
144 : #endif
145 676 : if (sector < 1 || sector > n) {
146 0 : AliWarning(Form("sector # %d out of range (0-%d)", sector-1, n-1));
147 0 : return kFALSE;
148 : }
149 338 : sector--;
150 : // Get track position
151 338 : mc->TrackPosition(v);
152 676 : AliFMDDebug(40, ("<2> Inside an active FMD volume FMD%d%c[%2d,%3d] %s",
153 : detector, ring, sector, strip, mc->CurrentVolPath()));
154 :
155 338 : return kTRUE;
156 338 : }
157 :
158 :
159 : //____________________________________________________________________
160 : Bool_t
161 : AliFMDv1::CheckHit(Int_t trackno, Int_t pdg, Float_t absQ,
162 : const TLorentzVector& p, Float_t edep) const
163 : {
164 : // Check that a hit is good
165 1014 : if (AliLog::GetDebugLevel("FMD", "AliFMD") < 5) return kFALSE;
166 0 : TVirtualMC* mc = TVirtualMC::GetMC();
167 0 : Double_t mass = mc->TrackMass();
168 0 : Double_t poverm = (mass == 0 ? 0 : p.P() / mass);
169 :
170 : // This `if' is to debug abnormal energy depositions. We trigger on
171 : // p/m approx larger than or equal to a MIP, and a large edep - more
172 : // than 1 keV - a MIP is 100 eV.
173 0 : if (!(edep > absQ * absQ && poverm > 1)) return kFALSE;
174 :
175 0 : TArrayI procs;
176 0 : mc->StepProcesses(procs);
177 0 : TString processes;
178 0 : for (Int_t ip = 0; ip < procs.fN; ip++) {
179 0 : if (ip != 0) processes.Append(",");
180 0 : processes.Append(TMCProcessName[procs.fArray[ip]]);
181 : }
182 0 : TDatabasePDG* pdgDB = TDatabasePDG::Instance();
183 0 : TParticlePDG* particleType = pdgDB->GetParticle(pdg);
184 0 : TString pname(particleType ? particleType->GetName() : "???");
185 0 : TString what;
186 0 : if (mc->IsTrackEntering()) what.Append("entering ");
187 0 : if (mc->IsTrackExiting()) what.Append("exiting ");
188 0 : if (mc->IsTrackInside()) what.Append("inside ");
189 0 : if (mc->IsTrackDisappeared()) what.Append("disappeared ");
190 0 : if (mc->IsTrackStop()) what.Append("stopped ");
191 0 : if (mc->IsNewTrack()) what.Append("new ");
192 0 : if (mc->IsTrackAlive()) what.Append("alive ");
193 0 : if (mc->IsTrackOut()) what.Append("out ");
194 :
195 0 : Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno);
196 0 : AliFMDDebug(15, ("Track # %5d deposits a lot of energy\n"
197 : " Volume: %s\n"
198 : " Momentum: (%7.4f,%7.4f,%7.4f)\n"
199 : " PDG: %d (%s)\n"
200 : " Edep: %-14.7f keV (mother %d)\n"
201 : " p/m: %-7.4f/%-7.4f = %-14.7f\n"
202 : " Processes: %s\n"
203 : " What: %s\n",
204 : trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
205 : pdg, pname.Data(), edep, mother, p.P(), mass,
206 : poverm, processes.Data(), what.Data()));
207 : return kTRUE;
208 338 : }
209 :
210 :
211 : //____________________________________________________________________
212 : void
213 : AliFMDv1::StepManager()
214 : {
215 : // Member function that is executed each time a hit is made in the
216 : // FMD. None-charged particles are ignored. Dead tracks are
217 : // ignored.
218 : //
219 : // The procedure is as follows:
220 : //
221 : // - IF NOT track is alive THEN RETURN ENDIF
222 : // - IF NOT particle is charged THEN RETURN ENDIF
223 : // - IF NOT volume name is "STRI" or "STRO" THEN RETURN ENDIF
224 : // - Get strip number (volume copy # minus 1)
225 : // - Get phi division number (mother volume copy #)
226 : // - Get module number (grand-mother volume copy #)
227 : // - section # = 2 * module # + phi division # - 1
228 : // - Get ring Id from volume name
229 : // - Get detector # from grand-grand-grand-mother volume name
230 : // - Get pointer to sub-detector object.
231 : // - Get track position
232 : // - IF track is entering volume AND track is inside real shape THEN
233 : // - Reset energy deposited
234 : // - Get track momentum
235 : // - Get particle ID #
236 : /// - ENDIF
237 : // - IF track is inside volume AND inside real shape THEN
238 : /// - Update energy deposited
239 : // - ENDIF
240 : // - IF track is inside real shape AND (track is leaving volume,
241 : // or it died, or it is stopped THEN
242 : // - Create a hit
243 : // - ENDIF
244 : //
245 21172 : TVirtualMC* mc = TVirtualMC::GetMC();
246 12544 : if (!mc->IsTrackAlive()) return;
247 8628 : Double_t absQ = TMath::Abs(mc->TrackCharge());
248 11089 : if (absQ <= 0) return;
249 :
250 6167 : Int_t copy;
251 6167 : Int_t vol = mc->CurrentVolID(copy);
252 6167 : AliFMDGeometry* fmd = AliFMDGeometry::Instance();
253 6167 : if (!fmd->IsActive(vol)) {
254 11658 : AliFMDDebug(50, ("Not an FMD volume %d '%s'",vol,mc->CurrentVolName()));
255 5829 : return;
256 : }
257 338 : TLorentzVector v;
258 338 : UShort_t detector;
259 338 : Char_t ring;
260 338 : UShort_t sector;
261 338 : UShort_t strip;
262 :
263 338 : if (fmd->IsDetailed()) {
264 676 : if (!VMC2FMD(copy, v, detector, ring, sector, strip)) return;
265 : } else {
266 0 : if (!VMC2FMD(v, detector, ring, sector, strip)) return;
267 : }
268 338 : TLorentzVector p;
269 338 : mc->TrackMomentum(p);
270 338 : Int_t trackno = gAlice->GetMCApp()->GetCurrentTrackNumber();
271 338 : Int_t pdg = mc->TrackPid();
272 676 : Double_t edep = mc->Edep() * 1000; // keV
273 676 : Bool_t isBad = CheckHit(trackno, pdg, absQ, p, edep);
274 :
275 : // Check that the track is actually within the active area
276 676 : Bool_t entering = mc->IsTrackEntering();
277 676 : Bool_t inside = mc->IsTrackInside();
278 1470 : Bool_t out = (mc->IsTrackExiting()|| mc->IsTrackDisappeared()||
279 228 : mc->IsTrackStop());
280 : // Reset the energy deposition for this track, and update some of
281 : // our parameters.
282 338 : if (entering) {
283 636 : AliFMDDebug(15, ("Track # %8d entering active FMD volume %s: "
284 : "Edep=%f (%f,%f,%f)", trackno, mc->CurrentVolPath(),
285 : edep, v.X(), v.Y(), v.Z()));
286 159 : fCurrentP = p;
287 159 : fCurrentV = v;
288 159 : fCurrentDeltaE = edep;
289 159 : fCurrentPdg = pdg; // mc->IdFromPDG(pdg);
290 159 : }
291 : // If the track is inside, then update the energy deposition
292 407 : if (inside && fCurrentDeltaE >= 0) {
293 69 : fCurrentDeltaE += edep;
294 276 : AliFMDDebug(15, ("Track # %8d inside active FMD volume %s: Edep=%f, "
295 : "Accumulated Edep=%f (%f,%f,%f)", trackno,
296 : mc->CurrentVolPath(), edep, fCurrentDeltaE,
297 : v.X(), v.Y(), v.Z()));
298 : }
299 : // The track exits the volume, or it disappeared in the volume, or
300 : // the track is stopped because it no longer fulfills the cuts
301 : // defined, then we create a hit.
302 338 : if (out) {
303 110 : if (fCurrentDeltaE >= 0) {
304 110 : fCurrentDeltaE += edep;
305 440 : AliFMDDebug(15, ("Track # %8d exiting active FMD volume %s: Edep=%g, "
306 : "Accumulated Edep=%g (%f,%f,%f)", trackno,
307 : mc->CurrentVolPath(), edep, fCurrentDeltaE,
308 : v.X(), v.Y(), v.Z()));
309 110 : TVector3 cur(v.Vect());
310 220 : cur -= fCurrentV.Vect();
311 110 : Double_t len = cur.Mag();
312 : AliFMDHit* h =
313 220 : AddHitByFields(trackno, detector, ring, sector, strip,
314 440 : fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(),
315 440 : fCurrentP.X(), fCurrentP.Y(), fCurrentP.Z(),
316 110 : fCurrentDeltaE, fCurrentPdg, fCurrentV.T(),
317 330 : len, mc->IsTrackDisappeared()||mc->IsTrackStop());
318 : // Add a copy
319 110 : if (isBad && fBad) {
320 0 : new ((*fBad)[fBad->GetEntries()]) AliFMDHit(*h);
321 : }
322 : // Check the geometry that we can get back the coordinates.
323 : #ifdef CHECK_TRANS
324 : Double_t x, y, z;
325 : fmd->Detector2XYZ(detector, ring, sector, strip, x, y ,z);
326 : AliFMDDebug(1, ("Hit at (%f,%f,%f), geometry says (%f,%f,%f)",
327 : fCurrentV.X(), fCurrentV.Y(), fCurrentV.Z(), x, y, z));
328 : #endif
329 110 : }
330 110 : fCurrentDeltaE = -1;
331 110 : }
332 17429 : }
333 : //___________________________________________________________________
334 : //
335 : // EOF
336 : //
|
