]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PHOS/AliPHOSv1.cxx
Starting a collection of QA/Comparison macros
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.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 /* 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 ClassImp(AliPHOSv1)
77
78 //____________________________________________________________________________
79 AliPHOSv1::AliPHOSv1() : fCPVDigits("AliPHOSCPVDigit",20)
80 {
81   //Def ctor.
82 }
83
84 //____________________________________________________________________________
85 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
86   AliPHOSv0(name,title), fCPVDigits("AliPHOSCPVDigit",20)
87 {
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   fHits= new TClonesArray("AliPHOSHit",1000) ;
103 //   fCPVDigits("AliPHOSCPVDigit",20);
104   gAlice->GetMCApp()->AddHitList(fHits) ; 
105
106   fNhits = 0 ;
107
108   fIshunt     =  2 ; // All hits are associated with primary particles
109 }
110
111 //____________________________________________________________________________
112 AliPHOSv1::~AliPHOSv1()
113 {
114   // dtor
115  if ( fHits) {
116     fHits->Delete() ; 
117     delete fHits ;
118     fHits = 0 ; 
119  }
120 }
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   AliPHOSGeometry * geom = GetGeometry() ; 
133
134   newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
135
136   for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
137     curHit = static_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
138     if(curHit->GetPrimary() != primary) break ; 
139            // We add hits with the same primary, while GEANT treats primaries succesively 
140     if( *curHit == *newHit ) {
141       *curHit + *newHit ;
142       deja = kTRUE ;
143     }
144   }
145          
146   if ( !deja ) {
147     new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
148     // get the block Id number
149     Int_t relid[4] ;
150     geom->AbsToRelNumbering(Id, relid) ;
151
152     fNhits++ ;
153   }
154   
155   delete newHit;
156 }
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 }
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   AliDetector::FinishEvent(); 
175 }
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   Int_t          relid[4] ;           // (box, layer, row, column) indices
182   Int_t          absid    ;           // absolute cell ID number
183   Float_t        xyzte[5]={-1000.,-1000.,-1000.,0.,0.}  ; // position wrt MRS, time and energy deposited
184   TLorentzVector pos      ;           // Lorentz vector of the track current position
185   Int_t          copy     ;
186
187   Int_t moduleNumber ;
188   
189   static Int_t idPCPQ = -1;
190   if (strstr(fTitle.Data(),"noCPV") == 0) 
191     idPCPQ = TVirtualMC::GetMC()->VolId("PCPQ");
192
193   if( TVirtualMC::GetMC()->CurrentVolID(copy) == idPCPQ &&
194      (TVirtualMC::GetMC()->IsTrackEntering() ) &&
195       TVirtualMC::GetMC()->TrackCharge() != 0) {      
196     
197     TVirtualMC::GetMC() -> TrackPosition(pos);
198
199     Float_t xyzm[3], xyzd[3] ;
200     Int_t i;
201     for (i=0; i<3; i++) xyzm[i] = pos[i];
202     TVirtualMC::GetMC() -> Gmtod (xyzm, xyzd, 1);    // transform coordinate from master to daughter system    
203
204
205     Float_t        xyd[3]={0,0,0}   ;   //local position of the entering
206     xyd[0]  = xyzd[0];
207     xyd[1]  =-xyzd[2];
208     xyd[2]  =-xyzd[1];
209     
210     // Current momentum of the hit's track in the local ref. system
211     TLorentzVector pmom     ;        //momentum of the particle initiated hit
212     TVirtualMC::GetMC() -> TrackMomentum(pmom);
213     Float_t pm[3], pd[3];
214     for (i=0; i<3; i++)  
215       pm[i]   = pmom[i];
216     
217     TVirtualMC::GetMC() -> Gmtod (pm, pd, 2);        // transform 3-momentum from master to daughter system
218     pmom[0] = pd[0];
219     pmom[1] =-pd[1];
220     pmom[2] =-pd[2];
221
222     // Digitize the current CPV hit:
223     
224     // 1. find pad response and    
225     TVirtualMC::GetMC()->CurrentVolOffID(3,moduleNumber);
226     moduleNumber--;
227     
228 //     TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0);   // array of digits for current hit
229     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     ndigits = fCPVDigits.GetEntriesFast();
239     for (idigit=0; idigit<ndigits-1; idigit++) {
240       AliPHOSCPVDigit  *cpvDigit1 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
241       Float_t x1 = cpvDigit1->GetXpad() ;
242       Float_t z1 = cpvDigit1->GetYpad() ;
243       for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
244         AliPHOSCPVDigit  *cpvDigit2 = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(jdigit));
245         Float_t x2 = cpvDigit2->GetXpad() ;
246         Float_t z2 = cpvDigit2->GetYpad() ;
247         if (x1==x2 && z1==z2) {
248           Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
249           cpvDigit2->SetQpad(qsumpad) ;
250           fCPVDigits.RemoveAt(idigit) ;
251         }
252       }
253     }
254     fCPVDigits.Compress() ;
255     
256     // 3. add digits to temporary hit list fTmpHits
257     
258     ndigits = fCPVDigits.GetEntriesFast();
259     for (idigit=0; idigit<ndigits; idigit++) {
260       AliPHOSCPVDigit  *cpvDigit = static_cast<AliPHOSCPVDigit*>(fCPVDigits.UncheckedAt(idigit));
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
271       xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
272       xyzte[4] = cpvDigit->GetQpad() ;                          // amplitude in a pad
273
274       Int_t primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
275       AddHit(fIshunt, primary, absid, xyzte);  
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();
281       }
282     }
283     fCPVDigits.Clear();
284   }
285
286  
287   static Int_t idPXTL = TVirtualMC::GetMC()->VolId("PXTL");  
288   if(TVirtualMC::GetMC()->CurrentVolID(copy) == idPXTL ) { //  We are inside a PBWO crystal
289
290     TVirtualMC::GetMC()->TrackPosition(pos) ;
291     xyzte[0] = pos[0] ;
292     xyzte[1] = pos[1] ;
293     xyzte[2] = pos[2] ;
294
295     Float_t lostenergy = TVirtualMC::GetMC()->Edep(); 
296     
297     //Put in the TreeK particle entering PHOS and all its parents
298     if ( TVirtualMC::GetMC()->IsTrackEntering() ){
299       Float_t xyzd[3] ;
300       TVirtualMC::GetMC() -> Gmtod (xyzte, xyzd, 1);    // transform coordinate from master to daughter system    
301       if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){   //Entered close to forward surface  
302         Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ; 
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         TVirtualMC::GetMC() -> Gmtod (vert, vertd, 1);    // transform coordinate from master to daughter system
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);
312           while ( parent != -1 ) {
313             part = gAlice->GetMCApp()->Particle(parent) ; 
314             part->SetBit(kKeepBit);
315             parent = part->GetFirstMother() ; 
316           }
317         }
318       }
319     }
320     if ( lostenergy != 0 ) {  // Track is inside the crystal and deposits some energy 
321       xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;     
322       
323       TVirtualMC::GetMC()->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
324       
325       Int_t strip ;
326       TVirtualMC::GetMC()->CurrentVolOffID(3, strip);
327       Int_t cell ;
328       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       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 :) 
337       absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() + 
338                     row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
339
340       
341       //Calculates the light yield, the number of photons produced in the
342       //crystal 
343       //There is no dependence of reponce on distance from energy deposition to APD
344       Float_t lightYield = gRandom->Poisson(AliPHOSSimParam::GetInstance()->GetLightFactor() * lostenergy) ;
345
346       //Calculates de energy deposited in the crystal  
347       xyzte[4] = AliPHOSSimParam::GetInstance()->GetAPDFactor() * lightYield ;
348       
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         }
364       }
365       else{
366         primary  =  gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() ); 
367       }
368       
369       // add current hit to the hit list
370       // Info("StepManager","%d %d", primary, tracknumber) ; 
371       AddHit(fIshunt, primary, absid, xyzte);
372         
373     } // there is deposited energy
374   } // we are inside a PHOS Xtal
375   
376 }
377
378 //____________________________________________________________________________
379 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
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
391   const Float_t kCelWr  = GetGeometry()->GetPadSizePhi()/2;  // Distance between wires (2 wires above 1 pad)
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
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();
410   Float_t eloss = kdEdx;
411
412   Float_t dZY   = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
413   Float_t dXY   = pX/pNorm * GetGeometry()->GetCPVGasThickness();
414   gRandom->Rannor(rnor1,rnor2);
415   eloss *= (1 + kDetR*rnor1) *
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;
419   Float_t zhit2 = zhit1 + dZY;
420   Float_t xhit2 = xhit1 + dXY;
421
422   Int_t   iwht1 = (Int_t) (xhit1 / kCelWr);           // wire (x) coordinate "in"
423   Int_t   iwht2 = (Int_t) (xhit2 / kCelWr);           // wire (x) coordinate "out"
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;
430     zxe[1][0] = (iwht1 + 0.5) * kCelWr;
431     zxe[2][0] =  eloss/2;
432     zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
433     zxe[1][1] = (iwht1 + 0.5) * kCelWr;
434     zxe[2][1] =  eloss/2;
435   }
436   else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
437     nIter = 3;
438     Int_t iwht3 = (iwht1 + iwht2) / 2;
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
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;
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 );
449     zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
450     zxe[1][0] =  xwht1;
451     zxe[2][0] =  eloss * egm1;
452     zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
453     zxe[1][1] =  xwht2;
454     zxe[2][1] =  eloss * egm2;
455     zxe[0][2] =  dXY*(xwht3-xwht1)/dXY + zhit1;
456     zxe[1][2] =  xwht3;
457     zxe[2][2] =  eloss * egm3;
458   }
459   else {                                   // incline 2-wire hit
460     nIter = 2;
461     Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
462     Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
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;
470     zxe[2][0] =  eloss * egm1;
471     zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
472     zxe[1][1] =  xwht2;
473     zxe[2][1] =  eloss * egm2;
474   }
475
476   // Finite size of ionization region
477
478   Int_t nCellZ  = GetGeometry()->GetNumberOfCPVPadsZ();
479   Int_t nCellX  = GetGeometry()->GetNumberOfCPVPadsPhi();
480   Int_t nz3     = (kNgamz+1)/2;
481   Int_t nx3     = (kNgamx+1)/2;
482   cpvDigits->Expand(nIter*kNgamx*kNgamz);
483   TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
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];
490     Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
491     Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
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;
498     for (Int_t iz=1; iz<=kNgamz; iz++) {
499       Int_t kzg = izcell + iz - nz3;
500       if (kzg<=0 || kzg>nCellZ) continue;
501       Float_t zg = (Float_t)(iz-nz3) - zc;
502       for (Int_t ix=1; ix<=kNgamx; ix++) {
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);
509         qpad += kNoise*rnor2;
510         if (qpad<0) continue;
511         
512         // Fill the array with pad response ID and amplitude
513         new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
514       }
515     }
516   }
517 }
518
519 //____________________________________________________________________________
520 Float_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
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();
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;
536 }
537
538 //____________________________________________________________________________
539 Double_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
550   const Double_t kA=1.0;
551   const Double_t kB=0.7;
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++) {
557     Double_t b1 = (2*i + 1) * kB;
558     cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
559   }
560   cumulPRF *= kA/(2*TMath::Pi());
561   return cumulPRF;
562 }
563