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 :
16 : /* $Id$ */
17 :
18 : /* History of cvs commits:
19 : *
20 : * $Log$
21 : * Revision 1.111 2007/07/24 09:41:19 morsch
22 : * AliStack included for kKeepBit.
23 : *
24 : * Revision 1.110 2007/03/10 08:58:52 kharlov
25 : * Protection for noCPV geometry
26 : *
27 : * Revision 1.109 2007/03/01 11:37:37 kharlov
28 : * Strip units changed from 8x1 to 8x2 (T.Pocheptsov)
29 : *
30 : * Revision 1.108 2007/02/02 09:40:50 alibrary
31 : * Includes required by ROOT head
32 : *
33 : * Revision 1.107 2007/02/01 10:34:47 hristov
34 : * Removing warnings on Solaris x86
35 : *
36 : * Revision 1.106 2006/11/14 17:11:15 hristov
37 : * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
38 : *
39 : * Revision 1.105 2006/09/13 07:31:01 kharlov
40 : * Effective C++ corrections (T.Pocheptsov)
41 : *
42 : * Revision 1.104 2005/05/28 14:19:05 schutz
43 : * Compilation warnings fixed by T.P.
44 : *
45 : */
46 :
47 : //_________________________________________________________________________
48 : // Implementation version v1 of PHOS Manager class
49 : //---
50 : //---
51 : // Layout EMC + CPV has name IHEP:
52 : // Produces hits for CPV, cumulated hits
53 : //---
54 : //---
55 : //*-- Author: Yves Schutz (SUBATECH)
56 :
57 :
58 : // --- ROOT system ---
59 : #include <TClonesArray.h>
60 : #include <TParticle.h>
61 : #include <TVirtualMC.h>
62 :
63 : // --- Standard library ---
64 :
65 :
66 : // --- AliRoot header files ---
67 : #include "AliPHOSCPVDigit.h"
68 : #include "AliPHOSGeometry.h"
69 : #include "AliPHOSHit.h"
70 : #include "AliPHOSv1.h"
71 : #include "AliRun.h"
72 : #include "AliMC.h"
73 : #include "AliStack.h"
74 : #include "AliPHOSSimParam.h"
75 :
76 20 : ClassImp(AliPHOSv1)
77 :
78 : //____________________________________________________________________________
79 24 : AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20)
80 60 : {
81 : //Def ctor.
82 24 : }
83 :
84 : //____________________________________________________________________________
85 : AliPHOSv1::AliPHOSv1(const char *name, const char *title):
86 2 : AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20)
87 5 : {
88 : //
89 : // We store hits :
90 : // - fHits (the "normal" one), which retains the hits associated with
91 : // the current primary particle being tracked
92 : // (this array is reset after each primary has been tracked).
93 : //
94 :
95 :
96 :
97 : // We do not want to save in TreeH the raw hits
98 : // But save the cumulated hits instead (need to create the branch myself)
99 : // It is put in the Digit Tree because the TreeH is filled after each primary
100 : // and the TreeD at the end of the event (branch is set in FinishEvent() ).
101 :
102 3 : fHits= new TClonesArray("AliPHOSHit",1000) ;
103 : // fCPVDigits("AliPHOSCPVDigit",20);
104 1 : gAlice->GetMCApp()->AddHitList(fHits) ;
105 :
106 1 : fNhits = 0 ;
107 :
108 1 : fIshunt = 2 ; // All hits are associated with primary particles
109 2 : }
110 :
111 : //____________________________________________________________________________
112 : AliPHOSv1::~AliPHOSv1()
113 78 : {
114 : // dtor
115 13 : if ( fHits) {
116 3 : fHits->Delete() ;
117 6 : delete fHits ;
118 3 : fHits = 0 ;
119 3 : }
120 39 : }
121 :
122 : //____________________________________________________________________________
123 : void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
124 : {
125 : // Add a hit to the hit list.
126 : // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
127 :
128 : Int_t hitCounter ;
129 : AliPHOSHit *newHit ;
130 : AliPHOSHit *curHit ;
131 : Bool_t deja = kFALSE ;
132 558606 : AliPHOSGeometry * geom = GetGeometry() ;
133 :
134 279303 : newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
135 :
136 37911098 : for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
137 12406769 : curHit = static_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
138 12406769 : if(curHit->GetPrimary() != primary) break ;
139 : // We add hits with the same primary, while GEANT treats primaries succesively
140 12406705 : if( *curHit == *newHit ) {
141 278938 : *curHit + *newHit ;
142 : deja = kTRUE ;
143 278938 : }
144 : }
145 :
146 279303 : if ( !deja ) {
147 365 : new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
148 : // get the block Id number
149 365 : Int_t relid[4] ;
150 365 : geom->AbsToRelNumbering(Id, relid) ;
151 :
152 365 : fNhits++ ;
153 365 : }
154 :
155 558606 : delete newHit;
156 279303 : }
157 :
158 : //____________________________________________________________________________
159 : void AliPHOSv1::FinishPrimary()
160 : {
161 : // called at the end of each track (primary) by AliRun
162 : // hits are reset for each new track
163 : // accumulate the total hit-multiplicity
164 :
165 224 : }
166 :
167 : //____________________________________________________________________________
168 : void AliPHOSv1::FinishEvent()
169 : {
170 : // called at the end of each event by AliRun
171 : // accumulate the hit-multiplicity and total energy per block
172 : // if the values have been updated check it
173 :
174 8 : AliDetector::FinishEvent();
175 4 : }
176 : //____________________________________________________________________________
177 : void AliPHOSv1::StepManager(void)
178 : {
179 : // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
180 :
181 911176 : Int_t relid[4] ; // (box, layer, row, column) indices
182 455588 : Int_t absid ; // absolute cell ID number
183 455588 : Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
184 455588 : TLorentzVector pos ; // Lorentz vector of the track current position
185 455588 : Int_t copy ;
186 :
187 455588 : Int_t moduleNumber ;
188 :
189 : static Int_t idPCPQ = -1;
190 455588 : if (fCreateCPV)
191 0 : idPCPQ = TVirtualMC::GetMC()->VolId("PCPQ");
192 :
193 1366764 : if( TVirtualMC::GetMC()->CurrentVolID(copy) == idPCPQ &&
194 0 : (TVirtualMC::GetMC()->IsTrackEntering() ) &&
195 0 : TVirtualMC::GetMC()->TrackCharge() != 0) {
196 :
197 0 : TVirtualMC::GetMC() -> TrackPosition(pos);
198 :
199 0 : Float_t xyzm[3], xyzd[3] ;
200 : Int_t i;
201 0 : for (i=0; i<3; i++) xyzm[i] = pos[i];
202 0 : TVirtualMC::GetMC() -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
203 :
204 :
205 0 : Float_t xyd[3]={0,0,0} ; //local position of the entering
206 0 : xyd[0] = xyzd[0];
207 0 : xyd[1] =-xyzd[2];
208 0 : xyd[2] =-xyzd[1];
209 :
210 : // Current momentum of the hit's track in the local ref. system
211 0 : TLorentzVector pmom ; //momentum of the particle initiated hit
212 0 : TVirtualMC::GetMC() -> TrackMomentum(pmom);
213 0 : Float_t pm[3], pd[3];
214 0 : for (i=0; i<3; i++)
215 0 : pm[i] = pmom[i];
216 :
217 0 : TVirtualMC::GetMC() -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
218 0 : pmom[0] = pd[0];
219 0 : pmom[1] =-pd[1];
220 0 : pmom[2] =-pd[2];
221 :
222 : // Digitize the current CPV hit:
223 :
224 : // 1. find pad response and
225 0 : TVirtualMC::GetMC()->CurrentVolOffID(3,moduleNumber);
226 0 : moduleNumber--;
227 :
228 : // TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
229 0 : CPVDigitize(pmom,xyd,&fCPVDigits);
230 :
231 : Float_t xmean = 0;
232 : Float_t zmean = 0;
233 : Float_t qsum = 0;
234 : Int_t idigit,ndigits;
235 :
236 : // 2. go through the current digit list and sum digits in pads
237 :
238 0 : ndigits = fCPVDigits.GetEntriesFast();
239 0 : for (idigit=0; idigit<ndigits-1; idigit++) {
240 0 : AliPHOSCPVDigit *cpvDigit1 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
241 0 : Float_t x1 = cpvDigit1->GetXpad() ;
242 0 : Float_t z1 = cpvDigit1->GetYpad() ;
243 0 : for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
244 0 : AliPHOSCPVDigit *cpvDigit2 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(jdigit));
245 0 : Float_t x2 = cpvDigit2->GetXpad() ;
246 0 : Float_t z2 = cpvDigit2->GetYpad() ;
247 0 : if (x1==x2 && z1==z2) {
248 0 : Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
249 0 : cpvDigit2->SetQpad(qsumpad) ;
250 0 : fCPVDigits.RemoveAt(idigit) ;
251 0 : }
252 : }
253 : }
254 0 : fCPVDigits.Compress() ;
255 :
256 : // 3. add digits to temporary hit list fTmpHits
257 :
258 0 : ndigits = fCPVDigits.GetEntriesFast();
259 0 : for (idigit=0; idigit<ndigits; idigit++) {
260 0 : AliPHOSCPVDigit *cpvDigit = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
261 0 : relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
262 0 : relid[1] =-1 ; // means CPV
263 0 : relid[2] = cpvDigit->GetXpad() ; // column number of a pad
264 0 : relid[3] = cpvDigit->GetYpad() ; // row number of a pad
265 :
266 : // get the absolute Id number
267 0 : GetGeometry()->RelToAbsNumbering(relid, absid) ;
268 :
269 : // add current digit to the temporary hit list
270 :
271 0 : xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
272 0 : xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
273 :
274 0 : Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
275 0 : AddHit(fIshunt, primary, absid, xyzte);
276 :
277 0 : if (cpvDigit->GetQpad() > 0.02) {
278 0 : xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
279 0 : zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
280 0 : qsum += cpvDigit->GetQpad();
281 0 : }
282 : }
283 0 : fCPVDigits.Clear();
284 0 : }
285 :
286 :
287 455592 : static Int_t idPXTL = TVirtualMC::GetMC()->VolId("PXTL");
288 1366764 : if(TVirtualMC::GetMC()->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
289 :
290 763524 : TVirtualMC::GetMC()->TrackPosition(pos) ;
291 763524 : xyzte[0] = pos[0] ;
292 763524 : xyzte[1] = pos[1] ;
293 763524 : xyzte[2] = pos[2] ;
294 :
295 1145286 : Float_t lostenergy = TVirtualMC::GetMC()->Edep();
296 :
297 : //Put in the TreeK particle entering PHOS and all its parents
298 1145286 : if ( TVirtualMC::GetMC()->IsTrackEntering() ){
299 77721 : Float_t xyzd[3] ;
300 155442 : TVirtualMC::GetMC() -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
301 233163 : if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
302 31 : Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
303 31 : TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
304 31 : Float_t vert[3],vertd[3] ;
305 31 : vert[0]=part->Vx() ;
306 31 : vert[1]=part->Vy() ;
307 31 : vert[2]=part->Vz() ;
308 62 : TVirtualMC::GetMC() -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
309 93 : if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
310 : //0.1 to get rid of numerical errors
311 16 : part->SetBit(kKeepBit);
312 172 : while ( parent != -1 ) {
313 70 : part = gAlice->GetMCApp()->Particle(parent) ;
314 70 : part->SetBit(kKeepBit);
315 70 : parent = part->GetFirstMother() ;
316 : }
317 : }
318 31 : }
319 77721 : }
320 381762 : if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
321 837909 : xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
322 :
323 558606 : TVirtualMC::GetMC()->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
324 :
325 279303 : Int_t strip ;
326 558606 : TVirtualMC::GetMC()->CurrentVolOffID(3, strip);
327 279303 : Int_t cell ;
328 558606 : TVirtualMC::GetMC()->CurrentVolOffID(2, cell);
329 :
330 : //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
331 : //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
332 : //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
333 837909 : Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
334 558606 : Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
335 :
336 : // Absid for 8x2-strips. Looks nice :)
337 1396515 : absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
338 1117212 : row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
339 :
340 : //Calculates the light yield, the number of photons produced in the
341 : //crystal
342 : //There is no dependence of reponce on distance from energy deposition to APD
343 837909 : Float_t lightYield = gRandom->Poisson(AliPHOSSimParam::GetInstance()->GetLightFactor() * lostenergy) ;
344 :
345 : //Calculates de energy deposited in the crystal
346 558606 : xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ;
347 :
348 : Int_t primary ;
349 558606 : if(fIshunt == 2){
350 558606 : primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
351 279303 : TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
352 4247876 : while ( !part->TestBit(kKeepBit) ) {
353 1844635 : primary = part->GetFirstMother() ;
354 3689270 : if(primary == -1){
355 1844635 : primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
356 0 : break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
357 : //surface of the crystal. In this case it may have no primary at all.
358 : //We can not easily separate this case from the case when this is part of the shower,
359 : //developed in the neighboring crystal.
360 : }
361 1844635 : part = gAlice->GetMCApp()->Particle(primary) ;
362 : }
363 279303 : }
364 : else{
365 0 : primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
366 : }
367 :
368 : // add current hit to the hit list
369 : // Info("StepManager","%d %d", primary, tracknumber) ;
370 279303 : AddHit(fIshunt, primary, absid, xyzte);
371 :
372 279303 : } // there is deposited energy
373 381762 : } // we are inside a PHOS Xtal
374 :
375 455588 : }
376 :
377 : //____________________________________________________________________________
378 : void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
379 : {
380 : // ------------------------------------------------------------------------
381 : // Digitize one CPV hit:
382 : // On input take exact 4-momentum p and position zxhit of the hit,
383 : // find the pad response around this hit and
384 : // put the amplitudes in the pads into array digits
385 : //
386 : // Author: Yuri Kharlov (after Serguei Sadovsky)
387 : // 2 October 2000
388 : // ------------------------------------------------------------------------
389 :
390 0 : const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
391 : const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
392 : const Float_t kdEdx = 4.0; // Average energy loss in CPV;
393 : const Int_t kNgamz = 5; // Ionization size in Z
394 : const Int_t kNgamx = 9; // Ionization size in Phi
395 : const Float_t kNoise = 0.03; // charge noise in one pad
396 :
397 0 : Float_t rnor1,rnor2;
398 :
399 : // Just a reminder on axes notation in the CPV module:
400 : // axis Z goes along the beam
401 : // axis X goes across the beam in the module plane
402 : // axis Y is a normal to the module plane showing from the IP
403 :
404 0 : Float_t hitX = zxhit[0];
405 0 : Float_t hitZ =-zxhit[1];
406 0 : Float_t pX = p.Px();
407 0 : Float_t pZ =-p.Pz();
408 0 : Float_t pNorm = p.Py();
409 : Float_t eloss = kdEdx;
410 :
411 0 : Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
412 0 : Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
413 0 : gRandom->Rannor(rnor1,rnor2);
414 0 : eloss *= (1 + kDetR*rnor1) *
415 0 : TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
416 0 : Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
417 0 : Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
418 0 : Float_t zhit2 = zhit1 + dZY;
419 0 : Float_t xhit2 = xhit1 + dXY;
420 :
421 0 : Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
422 0 : Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
423 :
424 : Int_t nIter;
425 0 : Float_t zxe[3][5];
426 0 : if (iwht1==iwht2) { // incline 1-wire hit
427 : nIter = 2;
428 0 : zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
429 0 : zxe[1][0] = (iwht1 + 0.5) * kCelWr;
430 0 : zxe[2][0] = eloss/2;
431 0 : zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
432 0 : zxe[1][1] = (iwht1 + 0.5) * kCelWr;
433 0 : zxe[2][1] = eloss/2;
434 0 : }
435 0 : else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
436 : nIter = 3;
437 0 : Int_t iwht3 = (iwht1 + iwht2) / 2;
438 0 : Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
439 0 : Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
440 0 : Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
441 0 : Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
442 0 : Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
443 0 : Float_t dxw1 = xhit1 - xwr13;
444 0 : Float_t dxw2 = xhit2 - xwr23;
445 0 : Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
446 0 : Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
447 0 : Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
448 0 : zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
449 0 : zxe[1][0] = xwht1;
450 0 : zxe[2][0] = eloss * egm1;
451 0 : zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
452 0 : zxe[1][1] = xwht2;
453 0 : zxe[2][1] = eloss * egm2;
454 0 : zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
455 0 : zxe[1][2] = xwht3;
456 0 : zxe[2][2] = eloss * egm3;
457 0 : }
458 : else { // incline 2-wire hit
459 : nIter = 2;
460 0 : Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
461 0 : Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
462 0 : Float_t xwr12 = (xwht1 + xwht2) / 2;
463 0 : Float_t dxw1 = xhit1 - xwr12;
464 0 : Float_t dxw2 = xhit2 - xwr12;
465 0 : Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
466 0 : Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
467 0 : zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
468 0 : zxe[1][0] = xwht1;
469 0 : zxe[2][0] = eloss * egm1;
470 0 : zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
471 0 : zxe[1][1] = xwht2;
472 0 : zxe[2][1] = eloss * egm2;
473 : }
474 :
475 : // Finite size of ionization region
476 :
477 0 : Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
478 0 : Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
479 : Int_t nz3 = (kNgamz+1)/2;
480 : Int_t nx3 = (kNgamx+1)/2;
481 0 : cpvDigits->Expand(nIter*kNgamx*kNgamz);
482 : TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
483 :
484 0 : for (Int_t iter=0; iter<nIter; iter++) {
485 :
486 0 : Float_t zhit = zxe[0][iter];
487 0 : Float_t xhit = zxe[1][iter];
488 0 : Float_t qhit = zxe[2][iter];
489 0 : Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
490 0 : Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
491 0 : if ( zcell<=0 || xcell<=0 ||
492 0 : zcell>=nCellZ || xcell>=nCellX) return;
493 0 : Int_t izcell = (Int_t) zcell;
494 0 : Int_t ixcell = (Int_t) xcell;
495 0 : Float_t zc = zcell - izcell - 0.5;
496 0 : Float_t xc = xcell - ixcell - 0.5;
497 0 : for (Int_t iz=1; iz<=kNgamz; iz++) {
498 0 : Int_t kzg = izcell + iz - nz3;
499 0 : if (kzg<=0 || kzg>nCellZ) continue;
500 0 : Float_t zg = (Float_t)(iz-nz3) - zc;
501 0 : for (Int_t ix=1; ix<=kNgamx; ix++) {
502 0 : Int_t kxg = ixcell + ix - nx3;
503 0 : if (kxg<=0 || kxg>nCellX) continue;
504 0 : Float_t xg = (Float_t)(ix-nx3) - xc;
505 :
506 : // Now calculate pad response
507 0 : Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
508 0 : qpad += kNoise*rnor2;
509 0 : if (qpad<0) continue;
510 :
511 : // Fill the array with pad response ID and amplitude
512 0 : new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
513 0 : }
514 0 : }
515 0 : }
516 0 : }
517 :
518 : //____________________________________________________________________________
519 : Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
520 : // ------------------------------------------------------------------------
521 : // Calculate the amplitude in one CPV pad using the
522 : // cumulative pad response function
523 : // Author: Yuri Kharlov (after Serguei Sadovski)
524 : // 3 October 2000
525 : // ------------------------------------------------------------------------
526 :
527 0 : Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
528 0 : Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
529 0 : Double_t z = zhit * GetGeometry()->GetPadSizeZ();
530 0 : Double_t x = xhit * GetGeometry()->GetPadSizePhi();
531 0 : Double_t amplitude = qhit *
532 0 : (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
533 0 : CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
534 0 : return (Float_t)amplitude;
535 : }
536 :
537 : //____________________________________________________________________________
538 : Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
539 : // ------------------------------------------------------------------------
540 : // Cumulative pad response function
541 : // It includes several terms from the CF decomposition in electrostatics
542 : // Note: this cumulative function is wrong since omits some terms
543 : // but the cell amplitude obtained with it is correct because
544 : // these omitting terms cancel
545 : // Author: Yuri Kharlov (after Serguei Sadovski)
546 : // 3 October 2000
547 : // ------------------------------------------------------------------------
548 :
549 : const Double_t kA=1.0;
550 : const Double_t kB=0.7;
551 :
552 0 : Double_t r2 = x*x + y*y;
553 0 : Double_t xy = x*y;
554 : Double_t cumulPRF = 0;
555 0 : for (Int_t i=0; i<=4; i++) {
556 0 : Double_t b1 = (2*i + 1) * kB;
557 0 : cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
558 : }
559 0 : cumulPRF *= kA/(2*TMath::Pi());
560 0 : return cumulPRF;
561 : }
562 :
|
