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