Removed wrong dependence on distance from energy deposition to APD
[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"
7587f5a5 74
75ClassImp(AliPHOSv1)
76
bea63bea 77//____________________________________________________________________________
02ab1add 78AliPHOSv1::AliPHOSv1():
43fbaae1 79 fLightYieldMean(0.),
80 fIntrinsicPINEfficiency(0.),
81 fLightYieldAttenuation(0.),
82 fRecalibrationFactor(0.),
83 fElectronsPerGeV(0.),
84 fAPDGain(0.),
85 fLightFactor(0.),
86 fAPDFactor(0.)
bea63bea 87{
43fbaae1 88 //Def ctor.
bea63bea 89}
90
7587f5a5 91//____________________________________________________________________________
92AliPHOSv1::AliPHOSv1(const char *name, const char *title):
43fbaae1 93 AliPHOSv0(name,title),
94 fLightYieldMean(0.),
95 fIntrinsicPINEfficiency(0.),
96 fLightYieldAttenuation(0.),
97 fRecalibrationFactor(0.),
98 fElectronsPerGeV(0.),
99 fAPDGain(0.),
100 fLightFactor(0.),
101 fAPDFactor(0.)
7587f5a5 102{
5f20d3fb 103 //
ed4205d8 104 // We store hits :
5f20d3fb 105 // - fHits (the "normal" one), which retains the hits associated with
106 // the current primary particle being tracked
107 // (this array is reset after each primary has been tracked).
108 //
fa412d9b 109
037cc66d 110
5f20d3fb 111
112 // We do not want to save in TreeH the raw hits
113 // But save the cumulated hits instead (need to create the branch myself)
114 // It is put in the Digit Tree because the TreeH is filled after each primary
7b326aac 115 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
5f20d3fb 116
ed4205d8 117 fHits= new TClonesArray("AliPHOSHit",1000) ;
5d12ce38 118 gAlice->GetMCApp()->AddHitList(fHits) ;
5f20d3fb 119
ed4205d8 120 fNhits = 0 ;
5f20d3fb 121
f6d1e5e1 122 fIshunt = 2 ; // All hits are associated with primary particles
7b326aac 123
9688c1dd 124 //Photoelectron statistics:
125 // The light yield is a poissonian distribution of the number of
126 // photons created in the PbWo4 crystal, calculated using following formula
127 // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
128 // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
129 // LightYieldMean is parameter calculated to be over 47000 photons per GeV
130 // APDEfficiency is 0.02655
131 // k_0 is 0.0045 from Valery Antonenko
132 // The number of electrons created in the APD is
133 // NumberOfElectrons = APDGain * LightYield
134 // The APD Gain is 300
135 fLightYieldMean = 47000;
136 fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
27f33ee5 137 fLightYieldAttenuation = 0.0045 ;
138 fRecalibrationFactor = 13.418/ fLightYieldMean ;
139 fElectronsPerGeV = 2.77e+8 ;
140 fAPDGain = 300. ;
141 fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
184569b0 142 fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
5f20d3fb 143}
144
7587f5a5 145//____________________________________________________________________________
bea63bea 146AliPHOSv1::~AliPHOSv1()
b2a60966 147{
bea63bea 148 // dtor
88cb7938 149 if ( fHits) {
ed4205d8 150 fHits->Delete() ;
151 delete fHits ;
152 fHits = 0 ;
184569b0 153 }
7587f5a5 154}
155
7587f5a5 156//____________________________________________________________________________
2af5445a 157void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
bea63bea 158{
159 // Add a hit to the hit list.
f6d1e5e1 160 // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
bea63bea 161
5f20d3fb 162 Int_t hitCounter ;
bea63bea 163 AliPHOSHit *newHit ;
5f20d3fb 164 AliPHOSHit *curHit ;
165 Bool_t deja = kFALSE ;
fa7cce36 166 AliPHOSGeometry * geom = GetGeometry() ;
bea63bea 167
2af5445a 168 newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
bea63bea 169
7854a24a 170 for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
29b077b5 171 curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
9688c1dd 172 if(curHit->GetPrimary() != primary) break ;
173 // We add hits with the same primary, while GEANT treats primaries succesively
ed4205d8 174 if( *curHit == *newHit ) {
f15a01eb 175 *curHit + *newHit ;
ed4205d8 176 deja = kTRUE ;
5f20d3fb 177 }
178 }
179
180 if ( !deja ) {
ed4205d8 181 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
7b326aac 182 // get the block Id number
9688c1dd 183 Int_t relid[4] ;
fa7cce36 184 geom->AbsToRelNumbering(Id, relid) ;
184569b0 185
ed4205d8 186 fNhits++ ;
5f20d3fb 187 }
184569b0 188
bea63bea 189 delete newHit;
bea63bea 190}
191
7b326aac 192//____________________________________________________________________________
193void AliPHOSv1::FinishPrimary()
194{
195 // called at the end of each track (primary) by AliRun
196 // hits are reset for each new track
197 // accumulate the total hit-multiplicity
7b326aac 198
199}
200
201//____________________________________________________________________________
202void AliPHOSv1::FinishEvent()
203{
204 // called at the end of each event by AliRun
205 // accumulate the hit-multiplicity and total energy per block
206 // if the values have been updated check it
88cb7938 207
88cb7938 208 AliDetector::FinishEvent();
7b326aac 209}
5f20d3fb 210//____________________________________________________________________________
7587f5a5 211void AliPHOSv1::StepManager(void)
212{
9688c1dd 213 // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
b2a60966 214
4f5bbbd4 215 Int_t relid[4] ; // (box, layer, row, column) indices
216 Int_t absid ; // absolute cell ID number
471193a8 217 Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
4f5bbbd4 218 TLorentzVector pos ; // Lorentz vector of the track current position
fa412d9b 219 Int_t copy ;
7587f5a5 220
9688c1dd 221 Int_t moduleNumber ;
222
39926a84 223 static Int_t idPCPQ = -1;
224 if (strstr(fTitle.Data(),"noCPV") == 0)
225 idPCPQ = gMC->VolId("PCPQ");
226
d6fb41ac 227 if( gMC->CurrentVolID(copy) == idPCPQ &&
9688c1dd 228 (gMC->IsTrackEntering() ) &&
229 gMC->TrackCharge() != 0) {
f6d1e5e1 230
9688c1dd 231 gMC -> TrackPosition(pos);
9225a7c4 232
9688c1dd 233 Float_t xyzm[3], xyzd[3] ;
234 Int_t i;
235 for (i=0; i<3; i++) xyzm[i] = pos[i];
9225a7c4 236 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
237
238
e3daf02c 239 Float_t xyd[3]={0,0,0} ; //local position of the entering
9688c1dd 240 xyd[0] = xyzd[0];
53e03a1e 241 xyd[1] =-xyzd[2];
242 xyd[2] =-xyzd[1];
f6d1e5e1 243
9688c1dd 244 // Current momentum of the hit's track in the local ref. system
245 TLorentzVector pmom ; //momentum of the particle initiated hit
246 gMC -> TrackMomentum(pmom);
247 Float_t pm[3], pd[3];
248 for (i=0; i<3; i++)
249 pm[i] = pmom[i];
f6d1e5e1 250
9688c1dd 251 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
252 pmom[0] = pd[0];
cf75bc19 253 pmom[1] =-pd[1];
254 pmom[2] =-pd[2];
f6d1e5e1 255
9688c1dd 256 // Digitize the current CPV hit:
257
258 // 1. find pad response and
259 gMC->CurrentVolOffID(3,moduleNumber);
260 moduleNumber--;
261
262 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
90cceaf6 263 CPVDigitize(pmom,xyd,cpvDigits);
fa412d9b 264
9688c1dd 265 Float_t xmean = 0;
266 Float_t zmean = 0;
267 Float_t qsum = 0;
268 Int_t idigit,ndigits;
269
270 // 2. go through the current digit list and sum digits in pads
271
272 ndigits = cpvDigits->GetEntriesFast();
273 for (idigit=0; idigit<ndigits-1; idigit++) {
29b077b5 274 AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
9688c1dd 275 Float_t x1 = cpvDigit1->GetXpad() ;
276 Float_t z1 = cpvDigit1->GetYpad() ;
277 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
29b077b5 278 AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
9688c1dd 279 Float_t x2 = cpvDigit2->GetXpad() ;
280 Float_t z2 = cpvDigit2->GetYpad() ;
281 if (x1==x2 && z1==z2) {
ad4aeaf4 282 Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
283 cpvDigit2->SetQpad(qsumpad) ;
9688c1dd 284 cpvDigits->RemoveAt(idigit) ;
fa412d9b 285 }
286 }
9688c1dd 287 }
288 cpvDigits->Compress() ;
289
290 // 3. add digits to temporary hit list fTmpHits
291
292 ndigits = cpvDigits->GetEntriesFast();
293 for (idigit=0; idigit<ndigits; idigit++) {
29b077b5 294 AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
9688c1dd 295 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
296 relid[1] =-1 ; // means CPV
297 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
298 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
299
300 // get the absolute Id number
301 GetGeometry()->RelToAbsNumbering(relid, absid) ;
302
303 // add current digit to the temporary hit list
304
471193a8 305 xyzte[3] = gMC->TrackTime() ;
306 xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
2af5445a 307
308 Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
309 AddHit(fIshunt, primary, absid, xyzte);
9688c1dd 310
311 if (cpvDigit->GetQpad() > 0.02) {
312 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
313 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
314 qsum += cpvDigit->GetQpad();
fa412d9b 315 }
fa412d9b 316 }
e534a69d 317 if (cpvDigits) {
318 cpvDigits->Delete();
319 delete cpvDigits;
320 cpvDigits=0;
321 }
9688c1dd 322 }
037cc66d 323
9688c1dd 324
d6fb41ac 325 static Int_t idPXTL = gMC->VolId("PXTL");
326 if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
9688c1dd 327
fa412d9b 328 gMC->TrackPosition(pos) ;
471193a8 329 xyzte[0] = pos[0] ;
330 xyzte[1] = pos[1] ;
331 xyzte[2] = pos[2] ;
597e6309 332
9688c1dd 333 Float_t lostenergy = gMC->Edep();
f6d1e5e1 334
335 //Put in the TreeK particle entering PHOS and all its parents
336 if ( gMC->IsTrackEntering() ){
337 Float_t xyzd[3] ;
471193a8 338 gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
97c3e101 339 if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
5d12ce38 340 Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
81d4c3d5 341 TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
342 Float_t vert[3],vertd[3] ;
343 vert[0]=part->Vx() ;
344 vert[1]=part->Vy() ;
345 vert[2]=part->Vz() ;
346 gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
2af5445a 347 if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
348 //0.1 to get rid of numerical errors
349 part->SetBit(kKeepBit);
81d4c3d5 350 while ( parent != -1 ) {
351 part = gAlice->GetMCApp()->Particle(parent) ;
81d4c3d5 352 part->SetBit(kKeepBit);
353 parent = part->GetFirstMother() ;
354 }
f6d1e5e1 355 }
356 }
357 }
9688c1dd 358 if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
471193a8 359 xyzte[3] = gMC->TrackTime() ;
f6d1e5e1 360
9688c1dd 361 gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
7b326aac 362
9688c1dd 363 Int_t strip ;
364 gMC->CurrentVolOffID(3, strip);
365 Int_t cell ;
366 gMC->CurrentVolOffID(2, cell);
9225a7c4 367
331c963c 368 //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
369 //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
370 //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
371 Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
372 Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
373
374 // Absid for 8x2-strips. Looks nice :)
9688c1dd 375 absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
331c963c 376 row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
377
9688c1dd 378
471193a8 379 //Calculates the light yield, the number of photons produced in the
9688c1dd 380 //crystal
9225a7c4 381 //There is no dependence of reponce on distance from energy deposition to APD
382 Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy) ;
471193a8 383
9688c1dd 384 //Calculates de energy deposited in the crystal
9225a7c4 385 xyzte[4] = fAPDFactor * lightYield ;
9688c1dd 386
2af5445a 387 Int_t primary ;
388 if(fIshunt == 2){
389 primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
390 TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
391 while ( !part->TestBit(kKeepBit) ) {
392 primary = part->GetFirstMother() ;
393 if(primary == -1){
394 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
395 break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
396 //surface of the crystal. In this case it may have no primary at all.
397 //We can not easily separate this case from the case when this is part of the shower,
398 //developed in the neighboring crystal.
399 }
400 part = gAlice->GetMCApp()->Particle(primary) ;
401 }
5a49626b 402 }
2af5445a 403 else
404 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
5a49626b 405
2af5445a 406
407
9688c1dd 408 // add current hit to the hit list
21cd0c07 409 // Info("StepManager","%d %d", primary, tracknumber) ;
2af5445a 410 AddHit(fIshunt, primary, absid, xyzte);
184569b0 411
fa412d9b 412 } // there is deposited energy
413 } // we are inside a PHOS Xtal
f6d1e5e1 414
fa412d9b 415}
416
417//____________________________________________________________________________
90cceaf6 418void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
fa412d9b 419{
420 // ------------------------------------------------------------------------
421 // Digitize one CPV hit:
422 // On input take exact 4-momentum p and position zxhit of the hit,
423 // find the pad response around this hit and
424 // put the amplitudes in the pads into array digits
425 //
426 // Author: Yuri Kharlov (after Serguei Sadovsky)
427 // 2 October 2000
428 // ------------------------------------------------------------------------
429
fa7cce36 430 const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
a3dfe79c 431 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
432 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
433 const Int_t kNgamz = 5; // Ionization size in Z
434 const Int_t kNgamx = 9; // Ionization size in Phi
435 const Float_t kNoise = 0.03; // charge noise in one pad
fa412d9b 436
437 Float_t rnor1,rnor2;
438
439 // Just a reminder on axes notation in the CPV module:
440 // axis Z goes along the beam
441 // axis X goes across the beam in the module plane
442 // axis Y is a normal to the module plane showing from the IP
443
444 Float_t hitX = zxhit[0];
445 Float_t hitZ =-zxhit[1];
446 Float_t pX = p.Px();
447 Float_t pZ =-p.Pz();
448 Float_t pNorm = p.Py();
a3dfe79c 449 Float_t eloss = kdEdx;
3d402178 450
21cd0c07 451//Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
7b326aac 452
fa7cce36 453 Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
454 Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
fa412d9b 455 gRandom->Rannor(rnor1,rnor2);
a3dfe79c 456 eloss *= (1 + kDetR*rnor1) *
fa7cce36 457 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
458 Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
459 Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
fa412d9b 460 Float_t zhit2 = zhit1 + dZY;
461 Float_t xhit2 = xhit1 + dXY;
462
a3dfe79c 463 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
464 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
fa412d9b 465
466 Int_t nIter;
467 Float_t zxe[3][5];
468 if (iwht1==iwht2) { // incline 1-wire hit
469 nIter = 2;
470 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
a3dfe79c 471 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
472 zxe[2][0] = eloss/2;
fa412d9b 473 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
a3dfe79c 474 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
475 zxe[2][1] = eloss/2;
fa412d9b 476 }
477 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
478 nIter = 3;
479 Int_t iwht3 = (iwht1 + iwht2) / 2;
a3dfe79c 480 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
481 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
482 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
fa412d9b 483 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
484 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
485 Float_t dxw1 = xhit1 - xwr13;
486 Float_t dxw2 = xhit2 - xwr23;
a3dfe79c 487 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
488 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
489 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
fa412d9b 490 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
491 zxe[1][0] = xwht1;
a3dfe79c 492 zxe[2][0] = eloss * egm1;
fa412d9b 493 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
494 zxe[1][1] = xwht2;
a3dfe79c 495 zxe[2][1] = eloss * egm2;
fa412d9b 496 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
497 zxe[1][2] = xwht3;
a3dfe79c 498 zxe[2][2] = eloss * egm3;
fa412d9b 499 }
500 else { // incline 2-wire hit
501 nIter = 2;
a3dfe79c 502 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
503 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
fa412d9b 504 Float_t xwr12 = (xwht1 + xwht2) / 2;
505 Float_t dxw1 = xhit1 - xwr12;
506 Float_t dxw2 = xhit2 - xwr12;
507 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
508 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
509 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
510 zxe[1][0] = xwht1;
a3dfe79c 511 zxe[2][0] = eloss * egm1;
fa412d9b 512 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
513 zxe[1][1] = xwht2;
a3dfe79c 514 zxe[2][1] = eloss * egm2;
fa412d9b 515 }
bea63bea 516
fa412d9b 517 // Finite size of ionization region
518
fa7cce36 519 Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
520 Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
a3dfe79c 521 Int_t nz3 = (kNgamz+1)/2;
522 Int_t nx3 = (kNgamx+1)/2;
523 cpvDigits->Expand(nIter*kNgamx*kNgamz);
29b077b5 524 TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
fa412d9b 525
526 for (Int_t iter=0; iter<nIter; iter++) {
527
528 Float_t zhit = zxe[0][iter];
529 Float_t xhit = zxe[1][iter];
530 Float_t qhit = zxe[2][iter];
fa7cce36 531 Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
532 Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
fa412d9b 533 if ( zcell<=0 || xcell<=0 ||
534 zcell>=nCellZ || xcell>=nCellX) return;
535 Int_t izcell = (Int_t) zcell;
536 Int_t ixcell = (Int_t) xcell;
537 Float_t zc = zcell - izcell - 0.5;
538 Float_t xc = xcell - ixcell - 0.5;
a3dfe79c 539 for (Int_t iz=1; iz<=kNgamz; iz++) {
fa412d9b 540 Int_t kzg = izcell + iz - nz3;
541 if (kzg<=0 || kzg>nCellZ) continue;
542 Float_t zg = (Float_t)(iz-nz3) - zc;
a3dfe79c 543 for (Int_t ix=1; ix<=kNgamx; ix++) {
fa412d9b 544 Int_t kxg = ixcell + ix - nx3;
545 if (kxg<=0 || kxg>nCellX) continue;
546 Float_t xg = (Float_t)(ix-nx3) - xc;
547
548 // Now calculate pad response
549 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
a3dfe79c 550 qpad += kNoise*rnor2;
fa412d9b 551 if (qpad<0) continue;
552
553 // Fill the array with pad response ID and amplitude
3d402178 554 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
fa412d9b 555 }
fa412d9b 556 }
fa412d9b 557 }
558}
559
560//____________________________________________________________________________
561Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
562 // ------------------------------------------------------------------------
563 // Calculate the amplitude in one CPV pad using the
564 // cumulative pad response function
565 // Author: Yuri Kharlov (after Serguei Sadovski)
566 // 3 October 2000
567 // ------------------------------------------------------------------------
568
fa7cce36 569 Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
570 Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
571 Double_t z = zhit * GetGeometry()->GetPadSizeZ();
572 Double_t x = xhit * GetGeometry()->GetPadSizePhi();
fa412d9b 573 Double_t amplitude = qhit *
574 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
575 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
576 return (Float_t)amplitude;
7587f5a5 577}
578
fa412d9b 579//____________________________________________________________________________
580Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
581 // ------------------------------------------------------------------------
582 // Cumulative pad response function
583 // It includes several terms from the CF decomposition in electrostatics
584 // Note: this cumulative function is wrong since omits some terms
585 // but the cell amplitude obtained with it is correct because
586 // these omitting terms cancel
587 // Author: Yuri Kharlov (after Serguei Sadovski)
588 // 3 October 2000
589 // ------------------------------------------------------------------------
590
a3dfe79c 591 const Double_t kA=1.0;
592 const Double_t kB=0.7;
fa412d9b 593
594 Double_t r2 = x*x + y*y;
595 Double_t xy = x*y;
596 Double_t cumulPRF = 0;
597 for (Int_t i=0; i<=4; i++) {
a3dfe79c 598 Double_t b1 = (2*i + 1) * kB;
fa412d9b 599 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
600 }
a3dfe79c 601 cumulPRF *= kA/(2*TMath::Pi());
fa412d9b 602 return cumulPRF;
603}
7eb9d12d 604