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